##Allen Schmaltz and Robert Schub
##Replication of Clare and Danilovic Paper


library(foreign)
library(xtable)
library(car)
library(robustbase)
library(mgcv)
library(MASS)
library(Zelig)
library(survey)
library(lme4)
library(VGAM)
library(mlogit)
library(sandwich)
library(gee)
library(Design)
library(lpSolve)
library(WhatIf)
library(lmtest)
library(zoo)

############################
## Load the data for replication of the
## four models in Table 1.
############################
library(foreign)
data <- read.dta(file="clare.dta")
head(data)
dim(data)

############################
## Table 1, Model 1
## The general pattern of the code for this first table
## is similar to that of the subsequent models in Table 1.
## The general model is run with Zelig, and then a separate
## approach is used to create the "clustered" standard errors (following the approach described 
## in "Cluster-robust standard errors using R" by Mahmood Arai, 2011).
## Note that many of these variables will be overwritten in running the subsequent
## models, but the returned value "ses" from coeftest() is saved with a unique name 
## (sesModel1, sesModel2, sesModel3, and sesModel4) for each model to be used to
## generate the table using xtable.
## Note that the variable name convention is different for the code replicating Table 2.
## Note that a number of the install.packages() calls have been commented out. You may need to 
## run those calls if you have not previously installed those packages.
############################

model1Vars <- c("cwinit", "dyad", "rivalry_count1", "weakmr_decline9", "weakmr2_decline9", "relcap", "jdem", "pol_rel", "pceyrs", "pceyrs2",  "pceyrs3")

dataNew <- data[, model1Vars]
head(dataNew)
nrow(dataNew)
dataNewNoNA <- na.omit(dataNew)
nrow(dataNewNoNA)
#12192
library(Zelig)
#install.packages("sandwich")
library(sandwich)
table1Model1 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = dataNewNoNA)
#(table1Model1$aic-18)/-2
#From Zelig documentation: Akaike's Information Criterion (minus twice the maximized log-likelihood plus twice the number of coefficients). Approximated Log-likelihood (comparable to the results from Stata):
#-3098.600


library(sandwich)
#install.packages("lmtest")
library(lmtest)

mat <- estfun(table1Model1)
mat <- na.omit(mat)
dim(mat)
N <- nrow(mat)
u <- apply(mat, 2, function(x) tapply(x, dataNewNoNA$dyad, sum))
u <- na.omit(u)
numberOfClusters <- length(unique(dataNewNoNA$dyad))
df <- (numberOfClusters / (numberOfClusters - 1))*((nrow(mat) - 1)/ (nrow(mat) - table1Model1$rank))
vcovCL <- df*sandwich(table1Model1, meat=crossprod(u)/N)
ses <- coeftest(table1Model1, vcovCL)
sesModel1 <- ses
summary(ses)



############################
## Table 1, Model 2
############################
table1Model2Vars <- c("cwinit", "dyad", "rivalry_count1", "weakmr_decline9", "nriv_weakdec9mr", "weakmr2_decline9", "relcap", "jdem", "pol_rel", "pceyrs", "pceyrs2",  "pceyrs3")

dataNew <- data[, table1Model2Vars]
head(dataNew)
nrow(dataNew)
dataNewNoNA <- na.omit(dataNew)
nrow(dataNewNoNA)
#12192
library(Zelig)
#install.packages("sandwich")
library(sandwich)
table1Model2 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = dataNewNoNA)
coef(table1Model2)


library(sandwich)
#install.packages("lmtest")
library(lmtest)

mat <- estfun(table1Model2)
mat <- na.omit(mat)
dim(mat)
N <- nrow(mat)
u <- apply(mat, 2, function(x) tapply(x, dataNewNoNA$dyad, sum))
u <- na.omit(u)
numberOfClusters <- length(unique(dataNewNoNA$dyad))
df <- (numberOfClusters / (numberOfClusters - 1))*((nrow(mat) - 1)/ (nrow(mat) - table1Model2$rank))
vcovCL <- df*sandwich(table1Model2, meat=crossprod(u)/N)
ses <- coeftest(table1Model2, vcovCL)
ses
sesModel2 <- ses

############################
## Table 1, Model 3
############################
table1Model3Vars <- c("force", "dyad", "rivalry_count1", "weakmr_decline9", "weakmr2_decline9", "relcap", "jdem", "pol_rel")

dataNew <- data[, table1Model3Vars]
head(dataNew)
nrow(dataNew)
dataNewNoNA <- na.omit(dataNew)
nrow(dataNewNoNA)
#946
library(Zelig)
#install.packages("sandwich")
library(sandwich)
table1Model3 <- zelig(force ~ rivalry_count1 + weakmr_decline9 + weakmr2_decline9 + relcap + jdem + pol_rel, model = "probit", data = dataNewNoNA)
coef(table1Model3)


library(sandwich)
#install.packages("lmtest")
library(lmtest)

mat <- estfun(table1Model3)
mat <- na.omit(mat)
dim(mat)
N <- nrow(mat)
u <- apply(mat, 2, function(x) tapply(x, dataNewNoNA$dyad, sum))
u <- na.omit(u)
numberOfClusters <- length(unique(dataNewNoNA$dyad))
df <- (numberOfClusters / (numberOfClusters - 1))*((nrow(mat) - 1)/ (nrow(mat) - table1Model3$rank))
vcovCL <- df*sandwich(table1Model3, meat=crossprod(u)/N)
ses <- coeftest(table1Model3, vcovCL)
ses
sesModel3 <- ses

############################
## Table 1, Model 4
############################
modelVars <- c("force", "dyad", "rivalry_count1", "weakmr_decline9", "nriv_weakdec9mr", "weakmr2_decline9", "relcap", "jdem", "pol_rel")

dataNew <- data[, modelVars]
head(dataNew)
nrow(dataNew)
dataNewNoNA <- na.omit(dataNew)
nrow(dataNewNoNA)
#946
library(Zelig)
#install.packages("sandwich")
library(sandwich)
modelOut <- zelig(force ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel, model = "probit", data = dataNewNoNA)
coef(modelOut)


library(sandwich)
#install.packages("lmtest")
library(lmtest)

mat <- estfun(modelOut)
mat <- na.omit(mat)
dim(mat)
N <- nrow(mat)
u <- apply(mat, 2, function(x) tapply(x, dataNewNoNA$dyad, sum))
u <- na.omit(u)
numberOfClusters <- length(unique(dataNewNoNA$dyad))
df <- (numberOfClusters / (numberOfClusters - 1))*((nrow(mat) - 1)/ (nrow(mat) - modelOut$rank))
vcovCL <- df*sandwich(modelOut, meat=crossprod(u)/N)
ses <- coeftest(modelOut, vcovCL)
ses
sesModel4 <- ses

############################
## Table 1, xtable
## Note that this will just produce the skeleton of the table
## with the main values. Further formatting was done in LaTeX. See the
## accompanying LaTeX code.
############################

#install.packages("xtable")
library(xtable)
xtable(summary(table1Model1), digits = 3)

table1a <- sesModel1[2:3, 1:2]
space <- c(0,0)
table1ab <- sesModel1[4:10, 1:2]
table1b <- sesModel1[1, 1:2]
table1 <- rbind(table1a, space, table1ab, table1b)
table1

table2a <- sesModel2[2:11, 1:2]
table2b <- sesModel2[1, 1:2]
table2 <- rbind(table2a, table2b)
table2

table3a <- sesModel3[2:3, 1:2]
space <- c(0,0)
table3ab <- sesModel3[4:7, 1:2]
table3b <- sesModel3[1, 1:2]
table3 <- rbind(table3a, space, table3ab, space, space, space, table3b)
table3

table4a <- sesModel4[2:8, 1:2]
space <- c(0,0)
table4b <- sesModel4[1, 1:2]
table4 <- rbind(table4a, space, space, space, table4b)
table4

fullTable1 <- cbind(table1, table2, table3, table4)
xtable(fullTable1, digits = 3)


##################################
### TABLE 2 ######################
##################################


##Load TABLE 2/Schub Data#########
cla<-read.dta(file="clare.dta")
head(cla)


##TABLE 2.1######################
head(cla)
cla2.1<-cla[,c(1,2,3,4,6,8,9,11,12,13,14,15,16,20)]

cla2.1.o<-na.omit(cla2.1)
dim(cla2.1.o)
head(cla2.1.o)

cla.pr.2.1<-zelig(cwrecip_corr ~ rivalry_count2 + weakmr2_decline9 + weakmr_decline9 + relcap + jdem + pol_rel, data=cla2.1.o, model="probit")
summary(cla.pr.2.1)


##Cluster Standard Errors 2.1#############
mat2.1<-estfun(cla.pr.2.1)
mat2.1<-na.omit(mat2.1)
dim(mat2.1)
N2.1<-nrow(mat2.1)
u2.1<-apply(mat2.1,2,function(x) tapply (x,cla2.1.o$dyad,sum))
u2.1<-na.omit(u2.1)
clust2.1<-length(unique(cla2.1.o$dyad))
a2.1<-clust2.1
b2.1<-nrow(mat2.1)
c2.1<-ncol(mat2.1)
df2.1<-((a2.1)/(a2.1-1))*((b2.1-1)/(b2.1-c2.1))
df2.1
vcovCL2.1<-df2.1*sandwich(cla.pr.2.1,meat=crossprod(u2.1)/N2.1)
(vcovCL2.1)
ses2.1<-coeftest(cla.pr.2.1,vcovCL2.1)
ses2.1

xtable(ses2.1[,1:2],digits=3)

##TABLE 2.2######################
head(cla)
cla2.2<-cla[,c(1,2,3,4,6,8,9,11,12,13,18,20)]

cla2.2.o<-na.omit(cla2.2)
dim(cla2.2.o)
head(cla2.2.o)

cla.pr.2.2<-zelig(cwrecip_corr ~ rivalry_count2 + weakmr2_decline9 + nriv2_weakdec9mr + weakmr_decline9 + relcap + jdem + pol_rel, data=cla2.2.o, model="probit")
summary(cla.pr.2.2)


##Cluster Standard Errors 2.2#############

mat2.2<-estfun(cla.pr.2.2)
mat2.2<-na.omit(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cla2.2.o$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cla2.2.o$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cla.pr.2.2,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cla.pr.2.2,vcovCL2.2)
ses2.2

xtable(ses2.2[,1:2],digits=3)

##TABLE 2.3######################
head(cla)
cla2.3<-cla[,c(1,2,3,4,6,8,9,11,12,13,18,19)]
head(cla2.3)
cla2.3init<-cla[cla[,5]==1,]
dim(cla2.3init)
cla2.3.o<-na.omit(cla2.3init)
dim(cla2.3.o)
head(cla2.3.o)

cla.pr.2.3<-zelig(force2 ~ rivalry_count2 + weakmr2_decline9 + weakmr_decline9 + relcap + jdem + pol_rel, data=cla2.3.o, model="probit")
summary(cla.pr.2.3)


##Cluster Standard Errors 2.3#############

mat2.3<-estfun(cla.pr.2.3)
mat2.3<-na.omit(mat2.3)
dim(mat2.3)
N2.3<-nrow(mat2.3)
u2.3<-apply(mat2.3,2,function(x) tapply (x,cla2.3.o$dyad,sum))
u2.3<-na.omit(u2.3)
clust2.3<-length(unique(cla2.3.o$dyad))
a2.3<-clust2.3
b2.3<-nrow(mat2.3)
c2.3<-ncol(mat2.3)
df2.3<-((a2.3)/(a2.3-1))*((b2.3-1)/(b2.3-c2.3))
df2.3
vcovCL2.3<-df2.3*sandwich(cla.pr.2.3,meat=crossprod(u2.3)/N2.3)
(vcovCL2.3)
ses2.3<-coeftest(cla.pr.2.3,vcovCL2.3)
ses2.3

xtable(ses2.3[,1:2],digits=3)

##TABLE 2.4################
head(cla)
cla2.4<-cla[,c(1,2,3,4,6,8,9,11,12,13,18,19)]
head(cla2.4)
cla2.4init<-cla[cla[,5]==1,]
dim(cla2.4init)
cla2.4.o<-na.omit(cla2.4init)
dim(cla2.4.o)
head(cla2.4.o)

cla.pr.2.4<-zelig(force2 ~ rivalry_count2 + weakmr2_decline9 + nriv2_weakdec9mr + weakmr_decline9 + relcap + jdem + pol_rel, data=cla2.4.o, model="probit")
summary(cla.pr.2.4)
names(summary(cla.pr.2.4))


##Cluster Standard Errors 2.4#############

mat2.4<-estfun(cla.pr.2.4)
mat2.4<-na.omit(mat2.4)
dim(mat2.4)
N2.4<-nrow(mat2.4)
u2.4<-apply(mat2.4,2,function(x) tapply (x,cla2.4.o$dyad,sum))
u2.4<-na.omit(u2.4)
clust2.4<-length(unique(cla2.4.o$dyad))
a2.4<-clust2.4
b2.4<-nrow(mat2.4)
c2.4<-ncol(mat2.4)
df2.4<-((a2.4)/(a2.4-1))*((b2.4-1)/(b2.4-c2.4))
df2.4
vcovCL2.4<-df2.4*sandwich(cla.pr.2.4,meat=crossprod(u2.4)/N2.4)
(vcovCL2.4)
ses2.4<-coeftest(cla.pr.2.4,vcovCL2.4)
ses2.4


xtable(ses2.4[,1:2],digits=3)


#### TABLE 2 AGGREGATED ####################
tab2.2a<-ses2.2[2:8,1:2]
tab2.2b<-ses2.2[1,1:2]
tab2.2<-rbind(tab2.2a,tab2.2b)
tab2.2

tab2.4a<-ses2.4[2:8,1:2]
tab2.4b<-ses2.4[1,1:2]
tab2.4<-rbind(tab2.4a,tab2.4b)
tab2.4

space<-c(0,0)
tab2.1a<-ses2.1[2:3,1:2]
tab2.1b<-ses2.1[4:7,1:2]
tab2.1c<-ses2.1[1,1:2]
tab2.1<-rbind(tab2.1a,space,tab2.1b,tab2.1c)
tab2.1

tab2.3a<-ses2.3[2:3,1:2]
tab2.3b<-ses2.3[4:7,1:2]
tab2.3c<-ses2.3[1,1:2]
tab2.3<-rbind(tab2.3a,space,tab2.3b,tab2.3c)
tab2.3

tab2.2<-cbind(tab2.1,tab2.2,tab2.3,tab2.4)
rownames(tab2.2)<-c("Number of Potential Opponents (Target)","Target's Weak Reputation","Number of Potential Opponents x Target's Weak Reputation","Challenger's Weak Reputation","Relative Capabilities","Joint Democracy","Politically Revelant","Constant")
xtable(tab2.2,digits=3,caption="Table 2 Reactive Reputation Building: Target's Reciprocation and Escalation (Probit Model)")

#################################################
### ROBUSTNESS TESTS TABLE 1 MODEL 4 ############
#################################################

###Construct Table1.4 Dataset##
X.A<-cbind(cla$rivalry_count1,cla$weakmr_decline9, cla$nriv_weakdec9mr,cla$weakmr2_decline9,cla$relcap,cla$jdem,cla$pol_rel)
dim(X.A)
ones<-rep(1,12626)
X.1<-as.matrix(cbind(ones,X.A))
y.1<-cla$force
xY<-as.matrix(na.omit(cbind(y.1,X.1)))
X<-xY[,-1]
y<-xY[,1]


##Create and optimize function for Table 1.4##
ll.prob1.4<-function(beta,y=y,X=X){
	apart<-pnorm(X%*%beta,log=TRUE)
	bpart<-pnorm(X%*%beta,log=TRUE,lower.tail=FALSE)
	logl<-sum(y*apart+(1-y)*bpart)
	return(logl)
	}
opt1.4<-optim(par=rep(0,8),fn=ll.prob1.4,y=y,X=X,method="BFGS",control=list(fnscale=-1,maxit=1000),hessian=TRUE)
opt1.4$par
ses <- sqrt(diag(solve(-opt1.4$hessian)))
table.dat <- cbind(opt1.4$par, ses)
table.dat

##Actual vs. Fitted Plot for Table 1.4##
fitted.Y <- pnorm(X%*%opt1.4$par)
Y<-y
bin.size <- .1 
mid <- seq(bin.size/2, 1-bin.size/2, .001)
accuracy <- c()
for (i in 1:length(mid)) {
  low <- mid[i] - bin.size/2
  high <- mid[i] + bin.size/2
  accuracy[i] <- mean(Y[fitted.Y < high & fitted.Y > low])
}

plot(mid, accuracy, 
     xlab="Mean Fitted Prob(Y=1)", ylab="Fraction of Y=1",
     ylim=c(0,1), xlim=c(0,1), main="Table 1 Model 4: Fitted Probabilities vs. Actual Outcomes")
lines(loess.smooth(mid, accuracy), col="blue", lwd=3)
abline(0, 1, lty=2)
rug(fitted.Y, ticksize = 0.01)   


pdf(file="InSamp.pdf")
plot(mid, accuracy,, 
     xlab="Mean Fitted Prob(Y=1)", ylab="Fraction of Y=1",
     ylim=c(0,1), xlim=c(0,1), main="Model 1.4: Fitted Probabilities vs. Actual Outcomes")
lines(loess.smooth(mid, accuracy), col="blue", lwd=3)
abline(0, 1, lty=2)
rug(fitted.Y, ticksize = 0.01)
dev.off()

cutoff <- c(seq(0, .1, by=.0001), seq(.1, 1, by=.005))

##ROC Plot for Table1.4##
correct0s <- c()
correct1s <- c()
for (i in 1:length(cutoff)) {
  predictions <- as.numeric(fitted.Y > cutoff[i])
  correct0s[i] <- mean(predictions[Y==0]==0)
  correct1s[i] <- mean(predictions[Y==1]==1)
}
plot(correct1s, correct0s, type="l", xlim=c(0,1), ylim=c(0,1), 
     main="Table 1 Model 4: ROC Curve",
	 xlab="Proportion of 1's correctly classified",
	 ylab="Proportion of 0's correctly classified")
abline(a=1, b=-1, lty=2)


pdf(file="InSampROC.pdf") 
plot(correct1s, correct0s, type="l", xlim=c(0,1), ylim=c(0,1), 
     main="Model 1.4: ROC Curve",
	 xlab="Proportion of 1's correctly classified",
	 ylab="Proportion of 0's correctly classified")
abline(a=1, b=-1, lty=2)
dev.off()


## Cross-validation of Fitted vs. Actual Table 1.4 Plot ##
##
set.seed(02138)
N <- length(Y) #how many observations do we have?
size <- floor(N/10) #generate our bins of observations
results <- list()

for (i in 1:10) {
  if(i==1) training <- (size+1):length(Y)
  if(i>1 && i <10) training <- (1:length(Y))[-c(((i-1)*size+1):((i)*size))] 
  if(i==10) training <- 1:(9*size)
  
  X.train <- X[training,]
  Y.train <- Y[training]
  X.test <- X[-training,]
  opt <- optim(par = rep(0,ncol(X)), fn = ll.prob1.4, y = Y.train, X=X.train, 
               method = "BFGS", control = list(fnscale = -1), hessian = TRUE) 
  fitted.Y <- pnorm(X.test%*%opt$par)			  
  results[[i]] <- fitted.Y	 
}

fitted.Y <- unlist(results)


bin.size <- .1 
mid <- seq(bin.size/2, 1-bin.size/2, .001)
accuracy <- c()
for (i in 1:length(mid)) {
  low <- mid[i] - bin.size/2
  high <- mid[i] + bin.size/2
  accuracy[i] <- mean(Y[fitted.Y < high & fitted.Y > low])
}
plot(mid, accuracy,, 
     xlab="Mean Fitted Prob(Y=1)", ylab="Fraction of Y=1",
     ylim=c(0,1), xlim=c(0,1), main="Cross Validation Model 1.4: Fitted")
lines(loess.smooth(mid, accuracy), col="blue", lwd=3)
abline(0, 1, lty=2)
rug(fitted.Y, ticksize = 0.01) 

pdf(file="OutSamp.pdf")
plot(mid, accuracy,, 
     xlab="Mean Fitted Prob(Y=1)", ylab="Fraction of Y=1",
     ylim=c(0,1), xlim=c(0,1), main="Cross Validation Model 1.4: Fitted")
lines(loess.smooth(mid, accuracy), col="blue", lwd=3)
abline(0, 1, lty=2)
rug(fitted.Y, ticksize = 0.01)
dev.off()
##FIT VS ACTUAL OUT OF SAMPLE IS HORRENDOUS###


## Cross Validation of ROC Table 1.4 Plot ##
cutoff <- c(seq(0, .1, by=.0001), seq(.1, 1, by=.005))
##Let's also do an ROC plot
correct0s <- c()
correct1s <- c()
for (i in 1:length(cutoff)) {
  predictions <- as.numeric(fitted.Y > cutoff[i])
  correct0s[i] <- mean(predictions[Y==0]==0)
  correct1s[i] <- mean(predictions[Y==1]==1)
}
plot(correct1s, correct0s, type="l", xlim=c(0,1), ylim=c(0,1), 
     main="Cross-Validation Model 1.4: ROC",
	 xlab="Proportion of 1's correctly classified",
	 ylab="Proportion of 0's correctly classified")
abline(a=1, b=-1, lty=2)


pdf(file="OutSampROC.pdf") 
plot(correct1s, correct0s, type="l", xlim=c(0,1), ylim=c(0,1), 
     main="Cross-Validation Model 1.4: ROC",
	 xlab="Proportion of 1's correctly classified",
	 ylab="Proportion of 0's correctly classified")
abline(a=1, b=-1, lty=2)
dev.off()
##ROC OUT OF SAMPLE IS TERRIBLE###



####################################
## CONVEX HULL TEST FOR MODEL 1.4 ##
####################################



##PREP WORK FOR THE CONVEX HULL TEST BELOW###
z.out<-zelig(force~rivalry_count1+weakmr_decline9+weakmr2_decline9+relcap+jdem+pol_rel,data=cl1.3,model="probit")
z.out$coefficients

cl1.3<-na.omit(cla[c("force","rivalry_count1","weakmr_decline9","weakmr2_decline9","relcap","jdem","pol_rel")])
dim(cl1.3)

cl1.4<-na.omit(cla[c("force","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","dyad")])
dim(cl1.3)

summary(z.out)

x.out<-setx(z.out)
s.out<-sim(z.out,x=x.out)
summary(s.out)
plot(s.out)

hist(cla1.3$weakmr_decline9)

x.low<-setx(z.out, weakmr_decline9=0)
x.high<-setx(z.out, weakmr_decline9=0.7)
s.out<- sim(z.out,x=x.low,x1=x.high)
summary(s.out)
plot(s.out)


##END PREP WORK###



## TEST OF CONVEX HULL TREATING WEAKNESS (weakmr_decline9) AS TREATMENT WHERE ZEROS ARE FLIPPED TO GREATER THAN ZERO AND THOSE GREATER THAN ZERO ARE FLIPPED TO ZERO TO CREATE THE COUNTERFACTUAL##

##HULL 1##
a<-(whatif(~ rivalry_count1+weakmr2_decline9+relcap+jdem+pol_rel,data=cl1.3[cl1.3$weakmr_decline9 == 0,],cfact=cl1.3[cl1.3$weakmr_decline9 >= 0.05,]))

##HULL 2--THE INVERSE CASE##
b<-(whatif(~ rivalry_count1+weakmr2_decline9+relcap+jdem+pol_rel,data=cl1.3[cl1.3$weakmr_decline9 >= 0.05,],cfact=cl1.3[cl1.3$weakmr_decline9 == 0,]))

##111/129 ARE IN HULL 1##
##467/817 ARE IN HULL 2##
578/946

##REMOVING OBSERVATIONS NOT IN CONVEX HULL 1##
cl1.30<-cl1.4[cl1.4[,3]==0,]
dim(cl1.30)
cl1.31<-cl1.4[cl1.4[,3]>=0.05,]
dim(cl1.31)
ah<-a$in.hull
head(cl1.30)
ah[ah==TRUE] <- 1
length(ah)
dim(cl1.31)
turn.off1<-cbind(ah,cl1.31)
turn.off1[1:36,1:2]
ah

##REMOVING OBSERVATIONS NOT IN CONVEX HULL 2##
bh<-b$in.hull
bh[bh==TRUE] <- 1
head(bh)
dim(cl1.30)
turn.off0<-cbind(bh,cl1.30)
turn.off0[1:20,1:2]

##CONSTRUCTING THE CONVEX HULL COMPLIANT DATASET##
sub1<-turn.off1[turn.off1[,1]==1,]
sub0<-turn.off0[turn.off0[,1]==1,]
sub1<-sub1[,-1]
sub0<-sub0[,-1]
cl.sub<-rbind(sub1,sub0)


dim(cl.sub)
head(cl.sub)



##check to make sure I kept the proper observations##
a.test<-(whatif(~ rivalry_count1+weakmr2_decline9+relcap+jdem+pol_rel,data=cl1.3[cl1.3$weakmr_decline9 == 0,],cfact=cl.sub[cl.sub$weakmr_decline9 >= 0.05,]))
summary(a.test)

b.test<-(whatif(~ rivalry_count1+weakmr2_decline9+relcap+jdem+pol_rel,data=cl1.3[cl1.3$weakmr_decline9 >= 0.05,],cfact=cl.sub[cl.sub$weakmr_decline9 == 0,]))
summary(b.test)
##end check##

###################################
###################################
#### MESS AROUND TIME #############
###################################
###################################

a<-(whatif(~ weakmr2_decline9+relcap+jdem+pol_rel,data=cl1.3[cl1.3$weakmr_decline9 == 0 & cl1.3$rivalry_count1 == 1,],cfact=cl1.3[cl1.3$weakmr_decline9 >= 0.05 & cl1.3$rivalry_count1 >=2 ,]))
summary(a)
##99/124 are in 

b<-(whatif(~ weakmr2_decline9+relcap+jdem+pol_rel,data=cl1.3[cl1.3$weakmr_decline9 >= 0.05 & cl1.3$rivalry_count1 >=2,],cfact=cl1.3[cl1.3$weakmr_decline9 == 0 & cl1.3$rivalry_count1==1,]))
summary(b)
##138/208 are in 

138+99
124+208
##TOTAL 237 of 332 ARE IN HULL##








##don't include rival count as a hull constraint##
a<-(whatif(~ weakmr2_decline9+relcap+jdem+pol_rel,data=cl1.3[cl1.3$weakmr_decline9 == 0,],cfact=cl1.3[cl1.3$weakmr_decline9 >= 0.05,]))
summary(a) 
#(123/129)
##HULL 2--THE INVERSE CASE##
b<-(whatif(~ weakmr2_decline9+relcap+jdem+pol_rel,data=cl1.3[cl1.3$weakmr_decline9 >= 0.05,],cfact=cl1.3[cl1.3$weakmr_decline9 == 0,]))
summary(b)
#(631/817)



#check balance#
library(Zelig)
library(MatchIt)
library(tcltk)
install.packages("cem")
library(cem)
head(cl1.3)

imb<-imbalance(group=cl1.3$weakmr_decline9,data=cl1.3, drop=c("force","weakmr_decline9","rivalry_count1"))


cl1.3$weak<-cl1.3$weakmr_decline9
cl1.3$weak[cl1.3$weak>0]<-1

cl1.3$rival<-cl1.3$rivalry_count1
cl1.3$rival[cl1.3$rival==1]<-0
cl1.3$rival[cl1.3$rival>1]<-1

imb<-imbalance(group=cl1.3$weak & cl1.3$rival,data=cl1.3, drop=c("force","weakmr_decline9","rivalry_count1","weak","rival"))
imb

hist(cl1.3$rival)




# TABLE 1 MODEL 2 #
head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)

test<-cl1.2[cl1.2$rivalry_count1 == 7 & cl1.2$weakmr_decline9>0.84,]
dim(test)
mean(test)

cl.ll<-cl1.2[cl1.2$rivalry_count1==1 & cl1.2$weakmr_decline9==0,]
cl.lh<-cl1.2[cl1.2$rivalry_count1==1 & cl1.2$weakmr_decline9>0,]
cl.hl<-cl1.2[cl1.2$rivalry_count1>1 & cl1.2$weakmr_decline9==0,]
cl.hh<-cl1.2[cl1.2$rivalry_count1>1 & cl1.2$weakmr_decline9>0,]
cl.non<-rbind(cl.ll,cl.lh,cl.hl)


head(cl.lh)
dim(cl.ll)

mean(cl.ll$cwinit)
mean(cl.lh$cwinit)
mean(cl.hl$cwinit)
mean(cl.hh$cwinit)

mean(cl.non$cwinit)

summary(cl1.2$rivalry_count1==1 & cl1.2$weakmr_decline9 == 0)
summary(cl1.2$rivalry_count1==1 & cl1.2$weakmr_decline9 > 0)
summary(cl1.2$rivalry_count1>1 & cl1.2$weakmr_decline9 == 0)
summary(cl1.2$rivalry_count1>1 & cl1.2$weakmr_decline9 > 0)
summary(cl1.2$rivalry_count1==1)
summary(cl1.2$rivalry_count1>1)
summary(cl1.2$weakmr_decline9==0)
summary(cl1.2$weakmr_decline9>0)

unique(cl1.2$weakmr_decline9)
summary(cl1.2$weakmr_decline9 > .7)

unique(cl1.2$rivalry_count1)
summary(cl1.2$rivalry_count1 > 3)

##################################
##DESIGNING DATASET OF INTEREST###
##################################

cl1.2$weak<-cl1.2$weakmr_decline9
cl1.2$weak[cl1.2$weakmr_decline9==0]<-0
cl1.2$weak[cl1.2$weakmr_decline9 > 0 & cl1.2$weakmr_decline9 < 0.7]<-1
cl1.2$weak[cl1.2$weakmr_decline9 > .7]<-2
hist(cl1.2$weak)

cl1.2$rival<-cl1.2$rivalry_count1
cl1.2$rival[cl1.2$rivalry_count1==1]<-0
cl1.2$rival[cl1.2$rivalry_count1 > 1 & cl1.2$rivalry_count1 < 4]<-1
cl1.2$rival[cl1.2$rivalry_count1 >3]<-2
hist(cl1.2$rival)

cl1.2$treat<-cl1.2$weak
cl1.2$treat<-NA
cl1.2$treat[cl1.2$weak==0 & cl1.2$rival==0]<-0
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==1]<-1
cl1.2$treat[cl1.2$weak==2 & cl1.2$rival==2]<-2

summary(cl1.2$treat)
dim(cl1.2)
mean(cl1.2$cwinit[cl1.2$weakmr_decline9>0 & cl1.2$rivalry_count1>1])
mean(cl1.2$cwinit[cl1.2$weakmr_decline9==0 & cl1.2$rivalry_count1==1])

cl1.2<-na.omit(cl1.2)
dim(cl1.2)
head(cl1.2)
hist(cl1.2$cwinit)


cl1.2no<-cl1.2[cl1.2$treat==0,]
cl1.2low<-cl1.2[cl1.2$treat==1,]
cl1.2high<-cl1.2[cl1.2$treat==2,]
dim(cl1.2high)
mean(cl1.2high$cwinit)

mean(cl1.2$cwinit[cl1.2$treat==0])
mean(cl1.2$cwinit[cl1.2$treat==1])
mean(cl1.2$cwinit[cl1.2$treat==2])
cl1.2x<-cl1.2[cl1.2$treat==0 | cl1.2$treat==1,]

table1Model2 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = cl1.2)
summary(table1Model2)


pre.imb<-imbalance(group=cl1.2$treat,data=cl1.2,drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","treat"))
pre.imb


##TRYING CEM##
auto.match<-matchit(formula=treat ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3,data=cl1.2x, method="cem")

head(cl1.2x)


#############################
## LETS TRY DICHOTOMOUS 1.2 #
#############################
cl1.2$weak<-cl1.2$weakmr_decline9
cl1.2$weak[cl1.2$weakmr_decline9==0]<-0
cl1.2$weak[cl1.2$weakmr_decline9 > 0]<-1
hist(cl1.2$weak)

cl1.2$rival<-cl1.2$rivalry_count1
cl1.2$rival[cl1.2$rivalry_count1==1]<-0
cl1.2$rival[cl1.2$rivalry_count1 > 1]<-1
hist(cl1.2$rival)

cl1.2$treat<-cl1.2$weak
cl1.2$treat<-NA
cl1.2$treat[cl1.2$weak==0 & cl1.2$rival==0]<-0
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==1]<-1

cl1.2<-na.omit(cl1.2)
dim(cl1.2)

user.match<-matchit(formula=treat ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3,data=cl1.2, method="cem")
summary(auto.match)
user.match
user.data<-match.data(user.match)
dim(user.data)


pre<-summary(user.match)$sum.all
post<-summary(user.match)$sum.match

pre.std<-(pre[,1]-pre[,2])/pre[,3]
post.std<-(post[,1]-post[,2])/post[,3]
table<-cbind(pre.std,post.std)
table


plot(x = pre.std, y = 1:nrow(pre),
     pch = 19, col = "burlywood2", xlab = "Standardized Imbalance: Propensity Score 1", axes = FALSE, 
     xlim = c(-3,3), ylab = "")
points(x = post.std, y = 1:nrow(post),
     pch = 19, col = "cadetblue2")
abline(v = 0, lty = "dashed")
axis(1)
axis(2, at = 1:8, labels = rownames(pre), las = 1, cex.axis = .8)
legend("topright", legend = c("Pre-match","Post-match"), pch = c(19,19),
   col = c("burlywood2","cadetblue2"), cex = .8)




head(cl1.2)
cem.match<-cem(treatment="treat",data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","treat"))

cem.match

z.out<-zelig(formula=cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=user.data, model="probit")
summary(cem1)

noncem<-zelig(formula=cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cla, model="probit")
summary(noncem)

noncem.lim<-zelig(formula=cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2, model="probit")
summary(noncem.lim)

unique(cla$dyad)

head(cla)
head(user.data)

mat2.2<-estfun(z.out)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,user.data$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(user.data$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(z.out,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(z.out,vcovCL2.2)
ses2.2

ses2.2[,1:2]

#######################################
## REAL ATTEMPT AT MATCHING FOR 1.2 ###
#######################################

#######################################
## create trichotomous treat ##########
## groups for each treatment ##########
#######################################


head(cla)

head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)

#Resetting rivalry count 1 and the interaction term
cl1.2$rivalry_count1_orig<-cl1.2$rivalry_count1
cl1.2$rivalry_count1<-cl1.2$rivalry_count1_orig-1

cl1.2$nriv_weakdec9mr_orig<-cl1.2$nriv_weakdec9mr
cl1.2$nriv_weakdec9mr<-cl1.2$rivalry_count1 * cl1.2$weakmr_decline9

#Drop original terms here
dim(cl1.2)
cl1.2<-cl1.2[,-14]
cl1.2<-cl1.2[,-13]
head(cl1.2)


head(cl1.2)

cl1.2$weak<-cl1.2$weakmr_decline9
cl1.2$weak[cl1.2$weakmr_decline9==0]<-0
cl1.2$weak[cl1.2$weakmr_decline9 > 0 & cl1.2$weakmr_decline9 < 0.7]<-1
cl1.2$weak[cl1.2$weakmr_decline9 > .7]<-2
hist(cl1.2$weak)

cl1.2$rival<-cl1.2$rivalry_count1
cl1.2$rival[cl1.2$rivalry_count1==0]<-0
cl1.2$rival[cl1.2$rivalry_count1 > 0 & cl1.2$rivalry_count1 < 3]<-1
cl1.2$rival[cl1.2$rivalry_count1 >2]<-2
hist(cl1.2$rival)



##summary of data ####################
summary(cl1.2$rival==2 & cl1.2$weak==2)
mean(cl1.2$cwinit[cl1.2$rival==0 & cl1.2$weak==0])


########################################
## IN ALL OF THE BELOW REGRESSIONS I HAVE TO PLUCK OUT VARIABLES THAT HAVE NOT VARIATION (WHICH EVER OF THE TREATMENTS IS FIXED AT ZERO PLUS THOSE THAT POST-MATCHING HAVE ONLY ONE VALUE LEFT) OTHERWISE THE SANDWICH FUNCTION CHOKES ####
#######################################



## TRICHOTOMOUS W/ UNCOARSENING IN FINAL MODEL ##



############################################################
################## MOVING RIVAL W/ WEAK AT 0 HAS NO EFFECT##
######## THIS IS THE RIGHT THING ########################
##let's try weighting but from subsetted w/ unmatched out##
cl1.2strong<-cl1.2[cl1.2$weak==0,]
dim(cl1.2strong)

cem.mat<-cem(treatment="rival", data=cl1.2strong, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.mat

head(cl1.2strong)
cl1.2strong$IN <- cem.mat$matched
weights<-cem.mat$w
cl1.2strong<-cbind(cl1.2strong,weights)
dim(cl1.2strong)

cl1.2strong$WA<-cl1.2strong$weights*cl1.2strong$cwinit
mean(cl1.2strong$WA[cl1.2strong$rival==0 & cl1.2strong$weights>0])
length(cl1.2lone$WA[cl1.2strong$weak==0 & cl1.2strong$weights>0])

cl1.2strongIN <- cl1.2strong[cl1.2strong$IN==TRUE,]
dim(cl1.2strongIN)

mean(cl1.2strongIN$WA[cl1.2strongIN$rival==2])
sd(cl1.2strongIN$WA[cl1.2strongIN$rival==2])
#sanity check--they are gaps in means but SEs are BIG!#

cem.strongW<-zelig(formula=cwinit ~  rivalry_count1 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2strongIN, model="probit", weights=cl1.2strongIN$weights)
summary(cem.strongW)
##getting error message but is generating something##


##sandwich chokes bc no variation in and weakmr2 so drop both from regression
mat2.2<-estfun(cem.strongW)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2strongIN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2strongIN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.strongW,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.strongW,vcovCL2.2)
ses2.2

##YIELDS ATE OF 0.0216 AND SE OF ..0338 FOR EFFECT OF RIVAL WHICH IS INSIGNIFICANT##
## OMG GETTING SAME EFFECT SIZE BUT DIFFERENT SEs (which mean insignificant! when subsetted) WHEN USING FULL SET WITH WEIGHTS VS. SUBSETTED SET DROPPING THE UNMATCHED OBS BUT USING THE SAME WEIGHTS (IE THE 0 WEIGHTS IN THE FULL SET AREN'T AKIN TO DROPPING THE OBSERVATIONS THEMSELVES)

##SO I THINK THE LARGER AND INSIGNIFICANT SEs ARE RIGHT BUT THEY NEED TO BE ADJUSTED FOR THE NON 1X1 MATCHING!!!





############################################################
################## MOVING WEAK W/ RIVAL AT 0 HAS NO EFFECT##
######## THIS IS THE RIGHT THING ########################
##let's try weighting but from subsetted w/ unmatched out##
cl1.2lone<-cl1.2[cl1.2$rival==0,]
dim(cl1.2lone)

cem.mat2<-cem(treatment="weak", data=cl1.2lone, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.mat2

head(cl1.2lone)
cl1.2lone$IN <- cem.mat2$matched
weights2<-cem.mat2$w
cl1.2lone<-cbind(cl1.2lone,weights2)
dim(cl1.2lone)

cl1.2lone$WA<-cl1.2lone$weights2*cl1.2lone$cwinit
mean(cl1.2lone$WA[cl1.2lone$weak==2 & cl1.2lone$weights2>0])
length(cl1.2lone$WA[cl1.2lone$weak==0 & cl1.2lone$weights2>0])

cl1.2loneIN <- cl1.2lone[cl1.2lone$IN==TRUE,]
dim(cl1.2loneIN)

mean(cl1.2loneIN$WA[cl1.2loneIN$weak==0])
sd(cl1.2loneIN$WA[cl1.2loneIN$weak==0])
#sanity check--they are gaps in means but SEs are BIG!#

cem.loneW<-zelig(formula=cwinit ~ weakmr_decline9 + relcap + jdem + pceyrs + pceyrs2 + pceyrs3, data=cl1.2loneIN, model="probit", weights=cl1.2loneIN$weights2)
summary(cem.loneW)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.loneW)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2loneIN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2loneIN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.loneW,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.loneW,vcovCL2.2)
ses2.2

##YIELDS ATE OF -0.289 AND SE OF .0395 FOR EFFECT OF WEAK WHICH IS HIGHLY SIGNIFICANT BUT IN WRONG DIRECTION##
## OMG GETTING SAME EFFECT SIZE BUT DIFFERENT SEs (.626 which mean insignificant! when subsetted) WHEN USING FULL SET WITH WEIGHTS VS. SUBSETTED SET DROPPING THE UNMATCHED OBS BUT USING THE SAME WEIGHTS (IE THE 0 WEIGHTS IN THE FULL SET AREN'T AKIN TO DROPPING THE OBSERVATIONS THEMSELVES)

##SO I THINK THE LARGER AND INSIGNIFICANT SEs ARE RIGHT BUT THEY NEED TO BE ADJUSTED FOR THE NON 1X1 MATCHING!!!


##################################################
##################################################
### CONTROL CASE SHOULD HAVE ALL 5 BOXES IN IT ###
##################################################
##################################################

##  WAIT ROB: if you are just going to uncoarsen the treatment brackets when running the parametric model why even bother breaking into 3 buckets rather than just dichotomizing then take the matched observations, uncoarsen and you get all the benefits of continuous treatment  #######

head(cl1.2)

##Generate treatment and control boxes (agg 5 boxes for control)
cl1.2$treat<-cl1.2$weak
cl1.2$treat[cl1.2$weak==0 | cl1.2$rival==0]<-0
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==1]<-1
cl1.2$treat[cl1.2$weak==2 & cl1.2$rival==1]<-2
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==2]<-3
cl1.2$treat[cl1.2$weak==2 & cl1.2$rival==2]<-4

table(cl1.2$treat)			
head(cl1.2)
dim(cl1.2)
mean(cl1.2$cwinit[cl1.2$treat==4])


##CHECKING OVERLAP ON COVARIATES BTW VARIOUS TREATMENT BOXES##
plot(density(x=cl1.2$pceyrs[cl1.2$treat==0]))
lines(density(x=cl1.2$pceyrs[cl1.2$treat==1]),col="red")
lines(density(x=cl1.2$pceyrs[cl1.2$treat==2]),col="green")
lines(density(x=cl1.2$pceyrs[cl1.2$treat==3]),col="blue")
lines(density(x=cl1.2$pceyrs[cl1.2$treat==4]),col="yellow")

##ONly thing found here is that the 0,0 box has many fewer weak rivals than the 2,2 box ##
plot(density(x=cl1.2$weakmr2_decline9[cl1.2$rival==0 & cl1.2$weak==0]))
lines(density(x=cl1.2$weakmr2_decline9[cl1.2$treat==4]),lty=2)

mean(cl1.2$weakmr2_decline9[cl1.2$treat==4])
mean(cl1.2$weakmr2_decline9[cl1.2$weak==0 & cl1.2$rival==0])



plot(density(x=cl1.2$jdem[cl1.2$rival==0 & cl1.2$weak==0]))
lines(density(x=cl1.2$jdem[cl1.2$treat==4]),lty=2)

mean(cl1.2$jdem[cl1.2$treat==4])
mean(cl1.2$jdem[cl1.2$weak==0 & cl1.2$rival==0])


plot(density(x=cl1.2$weakmr2_decline9[cl1.2$rival==0 & cl1.2$weak==0]))
lines(density(x=cl1.2$weakmr2_decline9[cl1.2$treat==4]),lty=2)



plot(density(x=cl1.2$pceyrs2[cl1.2$rival==0]),col="blue")
lines(density(x=cl1.2$pceyrs2[cl1.2$rival==1]),col="red")
lines(density(x=cl1.2$pceyrs2[cl1.2$rival==2]))

mean(cl1.2$pceyrs[cl1.2$treat==4])


##MATCH:: TRY ALL FIVE BOXES IN EACH STRATA##

head(cl1.2)

##try cutpoints to broaden##
weak2cut<-0.5
jdemcut<-0.35
relcapcut<-0.5
cuts<-list(weakmr2_decline9=weak2cut, jdem=jdemcut, relcap=relcapcut)
hist(cl1.2$relcap)

cem.mat3<-cem(treatment="treat", data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival"))
cem.mat3

##LEFT WITH 502C, 54T1, 33T2, 43T3, 44T4 w/ L1 of 0.636 with no specifcations##
502+54+33+43+44

bal<-imbalance(group=cl1.2$treat,data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","treat"))
bal

##UNBALANCED HAS L1 of 0.836##

head(cl1.2)
cl1.2$IN <- cem.mat3$matched
weights3<-cem.mat3$w
cl1.2<-cbind(cl1.2,weights3)
dim(cl1.2)

cl1.2$WA<-cl1.2$weights3*cl1.2$cwinit
mean(cl1.2$WA[cl1.2$treat==4 & cl1.2$weights3>0])
length(cl1.2$WA[cl1.2$treat==2 & cl1.2$weights3>0])

cl1.2IN <- cl1.2[cl1.2$IN==TRUE,]
dim(cl1.2IN)
head(cl1.2IN)

cl1.2IN[cl1.2IN$treat==2,]


mean(cl1.2IN$WA[cl1.2IN$treat==4])
sd(cl1.2IN$WA[cl1.2IN$treat==4])
#sanity check--they are gaps in means but SEs are BIG!#

head(cl1.2IN)

##pol_rel and weakmr2 are all the same in matched data so drop from regression#
cem.W<-zelig(formula=cwinit ~ weakmr_decline9 + rivalry_count1 + nriv_weakdec9mr + relcap + jdem + pceyrs + pceyrs2 + pceyrs3, data=cl1.2IN, model="probit", weights=cl1.2IN$weights3)
summary(cem.W)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.W)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.W,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.W,vcovCL2.2)
ses2.2


## SO HERE WE HAVE IT: INTERACTION TERM ISN'T SIGNIFICANT BUT LOWER ORDER TERM OF RIVALRY COUNT IS. INTERACTION TERM COEF IS -0.48 (SE 0.34) BUT SEs ARE WRONG BC OF MATCHING WITH REPLACEMENT ##


#################################################
### GENERATE PREDICTED PROBABILITES FOR EACH TREATMENT BOX THEN GRAPH THEM AS WE MOVE WEAKNESS AROUND ##

## Get median covariate values for whole dataset ##
head(cl1.2IN)
cov<-apply(X=cl1.2IN, MARGIN=2, FUN=median)
cov<-cov[c(6:7,9:11)]
cov
#create control covariates#
rival.con<-median(cl1.2IN$rivalry_count1[cl1.2IN$treat==0])
weak.con<-median(cl1.2IN$weakmr_decline9[cl1.2IN$treat==0])
int.con<-rival.con * weak.con
cov.con<-c(1,weak.con,rival.con,int.con,cov)
cov.con

rival.con<-median(cl1.2IN$rivalry_count1[cl1.2IN$treat==1])
weak.con<-median(cl1.2IN$weakmr_decline9[cl1.2IN$treat==1])
int.con<-rival.con * weak.con
cov.t1<-c(1,weak.con,rival.con,int.con,cov)
cov.t1

rival.con<-median(cl1.2IN$rivalry_count1[cl1.2IN$treat==2])
weak.con<-median(cl1.2IN$weakmr_decline9[cl1.2IN$treat==2])
int.con<-rival.con * weak.con
cov.t2<-c(1,weak.con,rival.con,int.con,cov)
cov.t2

rival.con<-median(cl1.2IN$rivalry_count1[cl1.2IN$treat==3])
weak.con<-median(cl1.2IN$weakmr_decline9[cl1.2IN$treat==3])
int.con<-rival.con * weak.con
cov.t3<-c(1,weak.con,rival.con,int.con,cov)
cov.t3

rival.con<-median(cl1.2IN$rivalry_count1[cl1.2IN$treat==4])
weak.con<-median(cl1.2IN$weakmr_decline9[cl1.2IN$treat==4])
int.con<-rival.con * weak.con
cov.t4<-c(1,weak.con,rival.con,int.con,cov)
cov.t4







names(cem.W)
beta<-cem.W$coefficients
beta
beta.draw<-mvrnorm(1000,mu=beta,Sigma=vcovCL2.2)

est.con<-pnorm(cov.con%*%t(beta.draw))
est.t1<-pnorm(cov.t1%*%t(beta.draw))
est.t2<-pnorm(cov.t2%*%t(beta.draw))
est.t3<-pnorm(cov.t3%*%t(beta.draw))
est.t4<-pnorm(cov.t4%*%t(beta.draw))


mean(est.con)
mean(est.t1)
mean(est.t2)
mean(est.t3)
mean(est.t4)

quantile(est.con,.025);quantile(est.con,.975)
quantile(est.t1,.025);quantile(est.t1,.975)
quantile(est.t2,.025);quantile(est.t2,.975)
quantile(est.t3,.025);quantile(est.t3,.975)
quantile(est.t3,.025);quantile(est.t4,.975)


##PLOT EACH ONE'S DISTRIBUTION##
plot(density(est.con))
lines(density(est.t4),lty=2)


lines(density(est.t2),lty=3)
lines(density(est.t3),lty=4)
lines(density(est.t4),lty=5)


##PLOT FD FOR EACH VS CONTROL##
t1<-est.t1-est.con
t2<-est.t2-est.con
t3<-est.t3-est.con
t4<-est.t4-est.con

plot(density(t1))
lines(density(t4),lty=2)
abline(v=mean(t1))
abline(v=mean(t4),lty=2)


###################################################
## FIGURES FROM PAPER: FLEXIBLE WITH PROPER DATA ##
###################################################
head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)

#Resetting rivalry count 1 and the interaction term
cl1.2$rivalry_count1_orig<-cl1.2$rivalry_count1
cl1.2$rivalry_count1<-cl1.2$rivalry_count1_orig-1

cl1.2$nriv_weakdec9mr_orig<-cl1.2$nriv_weakdec9mr
cl1.2$nriv_weakdec9mr<-cl1.2$rivalry_count1 * cl1.2$weakmr_decline9

#Drop original terms here
dim(cl1.2)
cl1.2<-cl1.2[,-14]
cl1.2<-cl1.2[,-13]
head(cl1.2)


head(cl1.2)

cl1.2$weak<-cl1.2$weakmr_decline9
cl1.2$weak[cl1.2$weakmr_decline9==0]<-0
cl1.2$weak[cl1.2$weakmr_decline9 > 0 & cl1.2$weakmr_decline9 < 0.7]<-1
cl1.2$weak[cl1.2$weakmr_decline9 > .7]<-2
hist(cl1.2$weak)

cl1.2$rival<-cl1.2$rivalry_count1
cl1.2$rival[cl1.2$rivalry_count1==0]<-0
cl1.2$rival[cl1.2$rivalry_count1 > 0 & cl1.2$rivalry_count1 < 3]<-1
cl1.2$rival[cl1.2$rivalry_count1 >2]<-2
hist(cl1.2$rival)


cl1.2$treat<-cl1.2$weak
cl1.2$treat[cl1.2$weak==0 | cl1.2$rival==0]<-0
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==1]<-1
cl1.2$treat[cl1.2$weak==2 & cl1.2$rival==1]<-2
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==2]<-3
cl1.2$treat[cl1.2$weak==2 & cl1.2$rival==2]<-4

cem.mat3<-cem(treatment="treat", data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival"))
cem.mat3

##LEFT WITH 502C, 54T1, 33T2, 43T3, 44T4##
502+54+33+43+44

head(cl1.2)
cl1.2$IN <- cem.mat3$matched
weights3<-cem.mat3$w
cl1.2<-cbind(cl1.2,weights3)
dim(cl1.2)

cl1.2$WA<-cl1.2$weights3*cl1.2$cwinit
mean(cl1.2$WA[cl1.2$treat==4 & cl1.2$weights3>0])
length(cl1.2$WA[cl1.2$treat==2 & cl1.2$weights3>0])

cl1.2IN <- cl1.2[cl1.2$IN==TRUE,]
dim(cl1.2IN)

cl1.2IN[cl1.2IN$treat==2,]


mean(cl1.2IN$WA[cl1.2IN$treat==4])
sd(cl1.2IN$WA[cl1.2IN$treat==4])
#sanity check--they are gaps in means but SEs are BIG!#

head(cl1.2IN)

##pol_rel and weakmr2 are all the same in matched data so drop from regression#
cem.W<-zelig(formula=cwinit ~ weakmr_decline9 + rivalry_count1 + nriv_weakdec9mr + relcap + jdem + pceyrs + pceyrs2 + pceyrs3, data=cl1.2IN, model="probit", weights=cl1.2IN$weights3)
summary(cem.W)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.W)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.W,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.W,vcovCL2.2)
ses2.2



names(cem.W)
beta<-cem.W$coefficients


beta.draw<-mvrnorm(1000,mu=beta,Sigma=vcovCL2.2)


##Create covariates plugging in proper rival number
dim(cl1.2IN)


beta
head(cl1.2IN)
cov<-apply(X=cl1.2IN, MARGIN=2, FUN=median)
cov<-cov[c(6:7,9:11)]
cov.con<-c(1,0,0,0,cov)

riv0<-cov.con
riv0[3]<-0

riv2<-cov.con
riv2[3]<-2

riv4<-cov.con
riv4[3]<-4

riv6<-cov.con
riv6[3]<-6

## Create range in weakness##
unique(cla$weakmr_decline9)
weak.range<-sort(unique(cla$weakmr_decline9), decreasing=TRUE)
length(weak.range)
weak.range

## Create matrix to hold it ##
sim0<-matrix(data=NA, ncol=length(weak.range),nrow=10000)
sim2<-matrix(data=NA, ncol=length(weak.range),nrow=10000)
sim4<-matrix(data=NA, ncol=length(weak.range),nrow=10000)
sim6<-matrix(data=NA, ncol=length(weak.range),nrow=10000)

length(riv0)

## Create for loop and allow BOTH weak and interaction to move ##
for (j in 1:length(weak.range)){
	riv0[2]<-weak.range[j]
	riv0[4]<-riv0[2] * riv0[3]
	sim0[,j]<-pnorm(riv0%*%t(beta.draw))
	}
for (j in 1:length(weak.range)){
	riv2[2]<-weak.range[j]
	riv2[4]<-riv2[2] * riv2[3]
	sim2[,j]<-pnorm(riv2%*%t(beta.draw))
	}
for (j in 1:length(weak.range)){
	riv4[2]<-weak.range[j]
	riv4[4]<-riv4[2] * riv4[3]
	sim4[,j]<-pnorm(riv4%*%t(beta.draw))
	}
for (j in 1:length(weak.range)){
	riv6[2]<-weak.range[j]
	riv6[4]<-riv6[2] * riv6[3]
	sim6[,j]<-pnorm(riv6%*%t(beta.draw))
	}

riv2[2]<-.9
riv2[4]<-riv2[2] * riv2[3]
riv2
((riv2%*%(beta.draw[1,])))

init0<-apply(sim0,2,mean)
init2<-apply(sim2,2,mean)
init4<-apply(sim4,2,mean)
init6<-apply(sim6,2,mean)

init0.low<-apply(sim0,2,quantile,0.025)
init0.high<-apply(sim0,2,quantile,0.975)
init2.low<-apply(sim2,2,quantile,0.025)
init2.high<-apply(sim2,2,quantile,0.975)	
init4.low<-apply(sim4,2,quantile,0.025)
init4.high<-apply(sim4,2,quantile,0.975)
init6.low<-apply(sim6,2,quantile,0.025)
init6.high<-apply(sim6,2,quantile,0.975)

range.recast<-c(1,2,3,4,5,6,7,8,9,10,11)
range.inv<-sort(weak.range,decreasing=TRUE)

##PLOT WITH PROPER MATCHED DATA ##
plot(range.recast,init0,type="l",lty=1, ylim=c(0,1))
lines(range.recast,init2,lty=2)
lines(range.recast,init4,lty=3)
lines(range.recast,init6,lty=4)
segments(x0=range.recast,x1=range.recast,y0=init6.low,y1=init6.high,lty=4)
legend("topleft",c("0 Rivals","2 Rivals","4 Rivals","6 Rivals"),lty=c(1,2,3,4))


init2
riv2[3]<-2
riv2[2]<-.9
riv2[4]<-riv2[2]*riv2[3]
length(riv2)
cof<-colMeans(beta.draw)
length(cof)
mean(pnorm(riv2%*%t(beta.draw)))





###################################################
## FIGURES FROM PAPER: FLEXIBLE WITH ALL DATA #####
###################################################
cla12 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = cla)
coef(cla12)

mat2.2<-estfun(cla12)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cla12,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cla12,vcovCL2.2)
ses2.2

beta<-ses2.2[,1]
beta

beta.draw<-mvrnorm(1000,mu=beta,Sigma=vcovCL2.2)
colMeans(beta.draw)
head(cl1.2)
cov<-apply(X=cl1.2, MARGIN=2, FUN=median)
cov<-cov[c(5:11)]
cov<-c(1,0,0,0,cov)
cov.con<-cov

##Create covariates plugging in proper rival number
dim(cl1.2)


cov.auth<-c(1,0,0,0,.9,.5,.36,1,16,684,42488)
cov.con<-cov.auth

riv0<-cov.con
riv0[2]<-0

riv2<-cov.con
riv2[2]<-2

riv4<-cov.con
riv4[2]<-4

riv6<-cov.con
riv6[2]<-6

## Create range in weakness##
unique(cla$weakmr_decline9)
weak.range<-sort(unique(cla$weakmr_decline9), decreasing=TRUE)
length(weak.range)
weak.range

## Create matrix to hold it ##
sim0<-matrix(data=NA, ncol=length(weak.range),nrow=10000)
sim2<-matrix(data=NA, ncol=length(weak.range),nrow=10000)
sim4<-matrix(data=NA, ncol=length(weak.range),nrow=10000)
sim6<-matrix(data=NA, ncol=length(weak.range),nrow=10000)

cbind(riv0,beta,riv0*beta)

length(riv0)
dim(beta.draw)

pnorm(t(riv0)%*%(beta.draw[1,]))

riv2
## Create for loop and allow BOTH weak and interaction to move ##
for (j in 1:length(weak.range)){
	riv0[3]<-weak.range[j]
	riv0[4]<-riv0[2] * riv0[3]
	sim0[,j]<-pnorm(riv0%*%t(beta.draw))
	}
	sim0
for (j in 1:length(weak.range)){
	riv2[3]<-weak.range[j]
	riv2[4]<-riv2[2] * riv2[3]
	sim2[,j]<-pnorm(riv2%*%t(beta.draw))
	}
for (j in 1:length(weak.range)){
	riv4[3]<-weak.range[j]
	riv4[4]<-riv4[2] * riv4[3]
	sim4[,j]<-pnorm(riv4%*%t(beta.draw))
	}
for (j in 1:length(weak.range)){
	riv6[3]<-weak.range[j]
	riv6[4]<-riv6[2] * riv6[3]
	sim6[,j]<-pnorm(riv6%*%t(beta.draw))
	}

init0<-apply(sim0,2,mean)
init2<-apply(sim2,2,mean)
init4<-apply(sim4,2,mean)
init6<-apply(sim6,2,mean)

init0
mean(sim0[,1])
quantile(sim0[,1],.025)
hist(sim0[,6])


init0.low<-apply(sim0,2,quantile,0.025)
init0.high<-apply(sim0,2,quantile,0.975)
init2.low<-apply(sim2,2,quantile,0.025)
init2.high<-apply(sim2,2,quantile,0.975)	
init4.low<-apply(sim4,2,quantile,0.025)
init4.high<-apply(sim4,2,quantile,0.975)
init6.low<-apply(sim6,2,quantile,0.025)
init6.high<-apply(sim6,2,quantile,0.975)	

init2

range.recast<-c(1,2,3,4,5,6,7,8,9,10,11)
range.inv<-sort(weak.range,decreasing=TRUE)

#### HERE IS THEIR PLOT WITH CI AROUND RIVALS =6 AND CHART MATCHES EXCEPT FOR WHOLE THING SEEMS HIGHER PROBABILITY THAN THE OTHERS FOUND ####
plot(range.recast,init0,type="l",lty=1,ylim=c(0.0,.2))
lines(range.recast,init2,lty=2)
lines(range.recast,init4,lty=3)
lines(range.recast,init6,lty=4)
segments(x0=range.recast,x1=range.recast,y0=init0.low,y1=init0.high,lty=1)
legend("topright",c("0 Rivals","2 Rivals","4 Rivals","6 Rivals"),lty=c(1,2,3,4))

colMeans(beta.draw)

#init6
#riv6[2]<-6
#riv6[3]<-.5
#riv6[4]<-riv6[2]*riv6[3]
#length(riv6)
#cof<-colMeans(beta.draw)
#length(cof)
#mean(pnorm(riv6%*%t(beta.draw)))


#######################################################
### MATCHED W/ 2x2 SETUP--3 CONTROL AND ONE TREATMENT##
#######################################################

head(cla)

head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)

#Resetting rivalry count 1 and the interaction term
cl1.2$rivalry_count1_orig<-cl1.2$rivalry_count1
cl1.2$rivalry_count1<-cl1.2$rivalry_count1_orig-1

cl1.2$nriv_weakdec9mr_orig<-cl1.2$nriv_weakdec9mr
cl1.2$nriv_weakdec9mr<-cl1.2$rivalry_count1 * cl1.2$weakmr_decline9

#Drop original terms here
dim(cl1.2)
cl1.2<-cl1.2[,-14]
cl1.2<-cl1.2[,-13]
head(cl1.2)


head(cl1.2)

cl1.2$treat<-cl1.2$weakmr_decline9
cl1.2$treat[cl1.2$weakmr_decline9==0 & cl1.2$rivalry_count1==0]<-0
cl1.2$treat[cl1.2$weakmr_decline9==0 & cl1.2$rivalry_count1>0]<-0
cl1.2$treat[cl1.2$weakmr_decline9>0 & cl1.2$rivalry_count1==0]<-0
cl1.2$treat[cl1.2$weakmr_decline9>0 & cl1.2$rivalry_count1>0]<-1
hist(cl1.2$treat)

head(cl1.2)
cem.mat3<-cem(treatment="treat", data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad"))
cem.mat3

##LEFT WITH 5105 C and 1222 T w/ L1 of 0.600 with no specifcations##
5105+1222

##Pre-balance##
pre.bal<-imbalance(group=cl1.2$treat, data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","treat"))
pre.bal

##PRE-BAL of L1=0.779##-->GET NICE CHART THIS WAY

head(cl1.2)
cl1.2$IN <- cem.mat3$matched
weights3<-cem.mat3$w
cl1.2<-cbind(cl1.2,weights3)
dim(cl1.2)

cl1.2$WA<-cl1.2$weights3*cl1.2$cwinit
mean(cl1.2$WA[cl1.2$treat==4 & cl1.2$weights3>0])
length(cl1.2$WA[cl1.2$treat==2 & cl1.2$weights3>0])

cl1.2IN <- cl1.2[cl1.2$IN==TRUE,]
dim(cl1.2IN)
head(cl1.2IN)

cl1.2IN[cl1.2IN$treat==2,]


mean(cl1.2IN$WA[cl1.2IN$treat==4])
sd(cl1.2IN$WA[cl1.2IN$treat==4])
#sanity check--they are gaps in means but SEs are BIG!#

head(cl1.2IN)

##pol_rel and weakmr2 are all the same in matched data so drop from regression#
cem.W<-zelig(formula=cwinit ~ weakmr_decline9 + rivalry_count1 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2IN, model="probit", weights=cl1.2IN$weights3)
summary(cem.W)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.W)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.W,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.W,vcovCL2.2)
ses2.2


## SO HERE WE HAVE IT: INTERACTION TERM ISN'T SIGNIFICANT IN DICHOTOMOUS RUN EITHER BUT INVERTED EFFECT IS SMALLER THAN IN TRICHOTOMOUS ATTEMPT: COEF -0.023 and SE 0.061: SO THIS IS A GENEROUS INTERPRETATION AND STILL NOT GETTING THERE ##











































############ WHAT IF I RUN THE MODEL WITH ALL UNITS MATCHING REJECTED###

head(cl1.2)
cl1.2$OUT <- cem.mat3$matched
weights3<-cem.mat3$w
cl1.2<-cbind(cl1.2,weights3)
dim(cl1.2)

cl1.2$WA<-cl1.2$weights3*cl1.2$cwinit
mean(cl1.2$WA[cl1.2$treat==4 & cl1.2$weights3>0])
length(cl1.2$WA[cl1.2$treat==2 & cl1.2$weights3>0])

cl1.2OUT <- cl1.2[cl1.2$OUT==FALSE,]
dim(cl1.2OUT)




mean(cl1.2OUT$cwinit[cl1.2OUT$treat==4])
sd(cl1.2IN$WA[cl1.2IN$treat==4])
#sanity check--they are gaps in means but SEs are BIG!#

head(cl1.2IN)

##pol_rel and weakmr2 are all the same in matched data so drop from #regression#
cem.out<-zelig(formula=cwinit ~ weakmr_decline9 + rivalry_count1 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2OUT, model="probit")
summary(cem.out)
##getting error message but is generating something##
plot(cem.out)



##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.out)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2OUT$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2OUT$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.out,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.out,vcovCL2.2)
ses2.2

#### IN UNMATCHED DATA, WE GET EXPECTED RESULT W/ SIGNIFICANCE ON THE INTERACTION TERM ####




#####################################################


## TRICHOTOMOUS W/ NO UNCOARSENING IN FINAL MODEL ##

head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)

#Resetting rivalry count 1 and the interaction term
cl1.2$rivalry_count1_orig<-cl1.2$rivalry_count1
cl1.2$rivalry_count1<-cl1.2$rivalry_count1_orig-1

cl1.2$nriv_weakdec9mr_orig<-cl1.2$nriv_weakdec9mr
cl1.2$nriv_weakdec9mr<-cl1.2$rivalry_count1 * cl1.2$weakmr_decline9

#Drop original terms here
dim(cl1.2)
cl1.2<-cl1.2[,-14]
cl1.2<-cl1.2[,-13]
head(cl1.2)


head(cl1.2)

cl1.2$weak<-cl1.2$weakmr_decline9
cl1.2$weak[cl1.2$weakmr_decline9==0]<-0
cl1.2$weak[cl1.2$weakmr_decline9 > 0 & cl1.2$weakmr_decline9 < 0.7]<-1
cl1.2$weak[cl1.2$weakmr_decline9 > .7]<-2
hist(cl1.2$weak)

cl1.2$rival<-cl1.2$rivalry_count1
cl1.2$rival[cl1.2$rivalry_count1==0]<-0
cl1.2$rival[cl1.2$rivalry_count1 > 0 & cl1.2$rivalry_count1 < 3]<-1
cl1.2$rival[cl1.2$rivalry_count1 >2]<-2
hist(cl1.2$rival)




############################################################
################## MOVING RIVAL W/ WEAK AT 0 HAS NO EFFECT##
######## THIS IS THE RIGHT THING ########################
##let's try weighting but from subsetted w/ unmatched out##







cl1.2strong<-cl1.2[cl1.2$weak==0,]
dim(cl1.2strong)

cem.mat<-cem(treatment="rival", data=cl1.2strong, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.mat

head(cl1.2strong)
cl1.2strong$IN <- cem.mat$matched
weights<-cem.mat$w
cl1.2strong<-cbind(cl1.2strong,weights)
dim(cl1.2strong)

cl1.2strong$WA<-cl1.2strong$weights*cl1.2strong$cwinit
mean(cl1.2strong$WA[cl1.2strong$rival==0 & cl1.2strong$weights>0])
length(cl1.2lone$WA[cl1.2strong$weak==0 & cl1.2strong$weights>0])

cl1.2strongIN <- cl1.2strong[cl1.2strong$IN==TRUE,]
dim(cl1.2strongIN)

mean(cl1.2strongIN$WA[cl1.2strongIN$rival==2])
sd(cl1.2strongIN$WA[cl1.2strongIN$rival==2])
#sanity check--they are gaps in means but SEs are BIG!#

cem.strongW<-zelig(formula=cwinit ~  rival + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2strongIN, model="probit", weights=cl1.2strongIN$weights)
summary(cem.strongW)
##getting error message but is generating something##


##sandwich chokes bc no variation in and weakmr2 so drop both from regression
mat2.2<-estfun(cem.strongW)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2strongIN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2strongIN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.strongW,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.strongW,vcovCL2.2)
ses2.2

##YIELDS ATE OF 0.0422 AND SE OF .0779 FOR EFFECT OF RIVAL WHICH IS INSIGNIFICANT##
## OMG GETTING SAME EFFECT SIZE BUT DIFFERENT SEs (which mean insignificant! when subsetted) WHEN USING FULL SET WITH WEIGHTS VS. SUBSETTED SET DROPPING THE UNMATCHED OBS BUT USING THE SAME WEIGHTS (IE THE 0 WEIGHTS IN THE FULL SET AREN'T AKIN TO DROPPING THE OBSERVATIONS THEMSELVES)

##SO I THINK THE LARGER AND INSIGNIFICANT SEs ARE RIGHT BUT THEY NEED TO BE ADJUSTED FOR THE NON 1X1 MATCHING!!!





############################################################
################## MOVING WEAK W/ RIVAL AT 0 HAS NO EFFECT##
######## THIS IS THE RIGHT THING ########################
##let's try weighting but from subsetted w/ unmatched out##
cl1.2lone<-cl1.2[cl1.2$rival==0,]
dim(cl1.2lone)

cem.mat2<-cem(treatment="weak", data=cl1.2lone, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.mat2

head(cl1.2lone)
cl1.2lone$IN <- cem.mat2$matched
weights2<-cem.mat2$w
cl1.2lone<-cbind(cl1.2lone,weights2)
dim(cl1.2lone)

cl1.2lone$WA<-cl1.2lone$weights2*cl1.2lone$cwinit
mean(cl1.2lone$WA[cl1.2lone$weak==2 & cl1.2lone$weights2>0])
length(cl1.2lone$WA[cl1.2lone$weak==0 & cl1.2lone$weights2>0])

cl1.2loneIN <- cl1.2lone[cl1.2lone$IN==TRUE,]
dim(cl1.2loneIN)

mean(cl1.2loneIN$WA[cl1.2loneIN$weak==0])
sd(cl1.2loneIN$WA[cl1.2loneIN$weak==0])
#sanity check--they are gaps in means but SEs are BIG!#

cem.loneW<-zelig(formula=cwinit ~ weak + relcap + jdem + pceyrs + pceyrs2 + pceyrs3, data=cl1.2loneIN, model="probit", weights=cl1.2loneIN$weights2)
summary(cem.loneW)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.loneW)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2loneIN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2loneIN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.loneW,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.loneW,vcovCL2.2)
ses2.2

##YIELDS ATE OF -0.106 AND SE OF .251 FOR EFFECT OF WEAK WHICH IS NOT SIGNIFICANT AND IN WRONG DIRECTION##


##SO I THINK THE LARGER AND INSIGNIFICANT SEs ARE RIGHT BUT THEY NEED TO BE ADJUSTED FOR THE NON 1X1 MATCHING!!!


##################################################
##################################################
### CONTROL CASE SHOULD HAVE ALL 5 BOXES IN IT ###
##################################################
##################################################


head(cl1.2)

##Generate treatment and control boxes (agg 5 boxes for control)
cl1.2$treat<-cl1.2$weak
cl1.2$treat[cl1.2$weak==0 | cl1.2$rival==0]<-0
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==1]<-1
cl1.2$treat[cl1.2$weak==2 & cl1.2$rival==1]<-2
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==2]<-3
cl1.2$treat[cl1.2$weak==2 & cl1.2$rival==2]<-4

table(cl1.2$treat)			
head(cl1.2)

##MATCH:: TRY ALL FIVE BOXES IN EACH STRATA##
cem.mat3<-cem(treatment="treat", data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival"))
cem.mat3

##LEFT WITH 502C, 54T1, 33T2, 43T3, 44T4##
502+54+33+43+44

head(cl1.2)
cl1.2$IN <- cem.mat3$matched
weights3<-cem.mat3$w
cl1.2<-cbind(cl1.2,weights3)
dim(cl1.2)

cl1.2$WA<-cl1.2$weights3*cl1.2$cwinit
mean(cl1.2$WA[cl1.2$treat==4 & cl1.2$weights3>0])
length(cl1.2$WA[cl1.2$treat==2 & cl1.2$weights3>0])

cl1.2IN <- cl1.2[cl1.2$IN==TRUE,]
dim(cl1.2IN)

mean(cl1.2IN$WA[cl1.2IN$treat==4])
sd(cl1.2IN$WA[cl1.2IN$treat==4])
#sanity check--they are gaps in means but SEs are BIG!#

head(cl1.2IN)

##pol_rel and weakmr2 are all the same in matched data so drop from regression#
cem.W<-zelig(formula=cwinit ~ weak + rival + treat + relcap + jdem + pceyrs + pceyrs2 + pceyrs3, data=cl1.2IN, model="probit", weights=cl1.2IN$weights3)
summary(cem.W)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.W)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.W,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.W,vcovCL2.2)
ses2.2


## SO HERE WE HAVE IT: IT'S ALL SIGNIFICANT WHEN WE DON'T UNCOARSEN ON THE TREATMENTS WITH INTERACTION COEF OF -0.44 (SE 0.22) BUT IN WRONG DIRECTION IT SEEMS?? ##

## CRAP, GET DIFFERENT RESULTS DEPENDING ON WHETHER TO UNCOARSEN THE TREATMENTS WHEN PUMPING MATCHED DATA THROUGH THE PARAMETRIC MODEL ##




###########################################################


## DICHOTOMOUS W/ UNCOARSENING IN FINAL MODEL ##

head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)

#Resetting rivalry count 1 and the interaction term
cl1.2$rivalry_count1_orig<-cl1.2$rivalry_count1
cl1.2$rivalry_count1<-cl1.2$rivalry_count1_orig-1

cl1.2$nriv_weakdec9mr_orig<-cl1.2$nriv_weakdec9mr
cl1.2$nriv_weakdec9mr<-cl1.2$rivalry_count1 * cl1.2$weakmr_decline9

#Drop original terms here
dim(cl1.2)
cl1.2<-cl1.2[,-14]
cl1.2<-cl1.2[,-13]
head(cl1.2)


head(cl1.2)

cl1.2$weak<-cl1.2$weakmr_decline9
cl1.2$weak[cl1.2$weakmr_decline9==0]<-0
cl1.2$weak[cl1.2$weakmr_decline9 > 0]<-1
hist(cl1.2$weak)

cl1.2$rival<-cl1.2$rivalry_count1
cl1.2$rival[cl1.2$rivalry_count1==0]<-0
cl1.2$rival[cl1.2$rivalry_count1 > 0]<-1
hist(cl1.2$rival)



############################################################
################## MOVING RIVAL W/ WEAK AT 0 HAS NO EFFECT##
######## THIS IS THE RIGHT THING ########################
##let's try weighting but from subsetted w/ unmatched out##
cl1.2strong<-cl1.2[cl1.2$weak==0,]
dim(cl1.2strong)

cem.mat<-cem(treatment="rival", data=cl1.2strong, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.mat

head(cl1.2strong)
cl1.2strong$IN <- cem.mat$matched
weights<-cem.mat$w
cl1.2strong<-cbind(cl1.2strong,weights)
dim(cl1.2strong)

cl1.2strong$WA<-cl1.2strong$weights*cl1.2strong$cwinit
mean(cl1.2strong$WA[cl1.2strong$rival==1 & cl1.2strong$weights>0])
length(cl1.2lone$WA[cl1.2strong$weak==0 & cl1.2strong$weights>0])

cl1.2strongIN <- cl1.2strong[cl1.2strong$IN==TRUE,]
dim(cl1.2strongIN)

mean(cl1.2strongIN$WA[cl1.2strongIN$rival==1])
sd(cl1.2strongIN$WA[cl1.2strongIN$rival==1])
#sanity check--they are gaps in means but SEs are BIG!#

cem.strongW<-zelig(formula=cwinit ~  rivalry_count1 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2strongIN, model="probit", weights=cl1.2strongIN$weights)
summary(cem.strongW)
##getting error message but is generating something##


##sandwich chokes bc no variation in and weakmr2 so drop both from regression
mat2.2<-estfun(cem.strongW)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2strongIN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2strongIN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.strongW,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.strongW,vcovCL2.2)
ses2.2

##YIELDS ATE OF -0.0018 AND SE OF .0236 FOR EFFECT OF RIVAL WHICH IS INSIGNIFICANT##

##SO I THINK THE LARGER AND INSIGNIFICANT SEs ARE RIGHT BUT THEY NEED TO BE ADJUSTED FOR THE NON 1X1 MATCHING!!!





############################################################
################## MOVING WEAK W/ RIVAL AT 0 HAS AN EFFECT##
######## THIS IS THE RIGHT THING ########################
##let's try weighting but from subsetted w/ unmatched out##
cl1.2lone<-cl1.2[cl1.2$rival==0,]
dim(cl1.2lone)

cem.mat2<-cem(treatment="weak", data=cl1.2lone, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.mat2

head(cl1.2lone)
cl1.2lone$IN <- cem.mat2$matched
weights2<-cem.mat2$w
cl1.2lone<-cbind(cl1.2lone,weights2)
dim(cl1.2lone)

cl1.2lone$WA<-cl1.2lone$weights2*cl1.2lone$cwinit
mean(cl1.2lone$WA[cl1.2lone$weak==1 & cl1.2lone$weights2>0])
length(cl1.2lone$WA[cl1.2lone$weak==0 & cl1.2lone$weights2>0])

cl1.2loneIN <- cl1.2lone[cl1.2lone$IN==TRUE,]
dim(cl1.2loneIN)

mean(cl1.2loneIN$WA[cl1.2loneIN$weak==0])
sd(cl1.2loneIN$WA[cl1.2loneIN$weak==0])
#sanity check--they are gaps in means but SEs are BIG!#

##no match for weakmr2 so drop that##
cem.loneW<-zelig(formula=cwinit ~ weakmr_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2loneIN, model="probit", weights=cl1.2loneIN$weights2)
summary(cem.loneW)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.loneW)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2loneIN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2loneIN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.loneW,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.loneW,vcovCL2.2)
ses2.2

##YIELDS ATE OF -0.909 AND SE OF .366 FOR EFFECT OF WEAK WHICH IS SIGNIFICANT BUT IN WRONG DIRECTION##


##SO I THINK THE LARGER AND INSIGNIFICANT SEs ARE RIGHT BUT THEY NEED TO BE ADJUSTED FOR THE NON 1X1 MATCHING!!!


##################################################
##################################################
### CONTROL CASE SHOULD HAVE ONLY 2 BOXES IN IT ##
##################################################
##################################################

##  WEAK HAS AN EFFECT  --> 2 BOXES IN CONTROL AND 2 TREATMENTS##

head(cl1.2)

##Generate treatment and control boxes (agg 5 boxes for control)
cl1.2$treat<-cl1.2$weak
cl1.2$treat[cl1.2$weak==0 & cl1.2$rival==0]<-0
cl1.2$treat[cl1.2$weak==0 & cl1.2$rival==1]<-0
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==0]<-1
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==1]<-2

table(cl1.2$treat)			
head(cl1.2)

##MATCH:: TRY ALL THREE BOXES IN EACH STRATA##
cem.mat3<-cem(treatment="treat", data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival"))
cem.mat3

##LEFT WITH 960C, 66T1, 190T2
960+66+190

head(cl1.2)
cl1.2$IN <- cem.mat3$matched
weights3<-cem.mat3$w
cl1.2<-cbind(cl1.2,weights3)
dim(cl1.2)

cl1.2$WA<-cl1.2$weights3*cl1.2$cwinit
mean(cl1.2$WA[cl1.2$treat==2 & cl1.2$weights3>0])
length(cl1.2$WA[cl1.2$treat==0 & cl1.2$weights3>0])

cl1.2IN <- cl1.2[cl1.2$IN==TRUE,]
dim(cl1.2IN)

mean(cl1.2IN$WA[cl1.2IN$treat==0])
sd(cl1.2IN$WA[cl1.2IN$treat==0])
#sanity check--they are gaps in means but SEs are BIG!#

head(cl1.2IN)

##pol_rel and weakmr2 are all the same in matched data so drop from regression#
cem.W<-zelig(formula=cwinit ~ weakmr_decline9 + rivalry_count1 + nriv_weakdec9mr + relcap + jdem + pceyrs + pceyrs2 + pceyrs3, data=cl1.2IN, model="probit", weights=cl1.2IN$weights3)
summary(cem.W)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.W)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.W,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.W,vcovCL2.2)
ses2.2


## SO HERE WE HAVE IT: INTERACTION TERM ISN'T SIGNIFICANT AND LOWER ORDER TERMS AREN'T EITHER. INTERACTION TERM COEF IS -0.0248 (SE 0.14) BUT SEs ARE WRONG BC OF MATCHING WITH REPLACEMENT ##


#####################################################



###########################################################


## DICHOTOMOUS W/ NO COARSENING IN FINAL MODEL ##

head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)

#Resetting rivalry count 1 and the interaction term
cl1.2$rivalry_count1_orig<-cl1.2$rivalry_count1
cl1.2$rivalry_count1<-cl1.2$rivalry_count1_orig-1

cl1.2$nriv_weakdec9mr_orig<-cl1.2$nriv_weakdec9mr
cl1.2$nriv_weakdec9mr<-cl1.2$rivalry_count1 * cl1.2$weakmr_decline9

#Drop original terms here
dim(cl1.2)
cl1.2<-cl1.2[,-14]
cl1.2<-cl1.2[,-13]
head(cl1.2)


head(cl1.2)

cl1.2$weak<-cl1.2$weakmr_decline9
cl1.2$weak[cl1.2$weakmr_decline9==0]<-0
cl1.2$weak[cl1.2$weakmr_decline9 > 0]<-1
hist(cl1.2$weak)

cl1.2$rival<-cl1.2$rivalry_count1
cl1.2$rival[cl1.2$rivalry_count1==0]<-0
cl1.2$rival[cl1.2$rivalry_count1 > 0]<-1
hist(cl1.2$rival)



############################################################
################## MOVING RIVAL W/ WEAK AT 0 HAS NO EFFECT##
######## THIS IS THE RIGHT THING ########################
##let's try weighting but from subsetted w/ unmatched out##
cl1.2strong<-cl1.2[cl1.2$weak==0,]
dim(cl1.2strong)

cem.mat<-cem(treatment="rival", data=cl1.2strong, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.mat

head(cl1.2strong)
cl1.2strong$IN <- cem.mat$matched
weights<-cem.mat$w
cl1.2strong<-cbind(cl1.2strong,weights)
dim(cl1.2strong)

cl1.2strong$WA<-cl1.2strong$weights*cl1.2strong$cwinit
mean(cl1.2strong$WA[cl1.2strong$rival==1 & cl1.2strong$weights>0])
length(cl1.2lone$WA[cl1.2strong$weak==0 & cl1.2strong$weights>0])

cl1.2strongIN <- cl1.2strong[cl1.2strong$IN==TRUE,]
dim(cl1.2strongIN)

mean(cl1.2strongIN$WA[cl1.2strongIN$rival==1])
sd(cl1.2strongIN$WA[cl1.2strongIN$rival==1])
#sanity check--they are gaps in means but SEs are BIG!#

cem.strongW<-zelig(formula=cwinit ~  rival + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2strongIN, model="probit", weights=cl1.2strongIN$weights)
summary(cem.strongW)
##getting error message but is generating something##


##sandwich chokes bc no variation in and weakmr2 so drop both from regression
mat2.2<-estfun(cem.strongW)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2strongIN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2strongIN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.strongW,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.strongW,vcovCL2.2)
ses2.2

##YIELDS ATE OF -0.0146 AND SE OF .0851 FOR EFFECT OF RIVAL WHICH IS INSIGNIFICANT##

##SO I THINK THE LARGER AND INSIGNIFICANT SEs ARE RIGHT BUT THEY NEED TO BE ADJUSTED FOR THE NON 1X1 MATCHING!!!





############################################################
################## MOVING WEAK W/ RIVAL AT 0 HAS AN EFFECT##
######## THIS IS THE RIGHT THING ########################
##let's try weighting but from subsetted w/ unmatched out##
cl1.2lone<-cl1.2[cl1.2$rival==0,]
dim(cl1.2lone)

cem.mat2<-cem(treatment="weak", data=cl1.2lone, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.mat2

head(cl1.2lone)
cl1.2lone$IN <- cem.mat2$matched
weights2<-cem.mat2$w
cl1.2lone<-cbind(cl1.2lone,weights2)
dim(cl1.2lone)

cl1.2lone$WA<-cl1.2lone$weights2*cl1.2lone$cwinit
mean(cl1.2lone$WA[cl1.2lone$weak==1 & cl1.2lone$weights2>0])
length(cl1.2lone$WA[cl1.2lone$weak==0 & cl1.2lone$weights2>0])

cl1.2loneIN <- cl1.2lone[cl1.2lone$IN==TRUE,]
dim(cl1.2loneIN)

mean(cl1.2loneIN$WA[cl1.2loneIN$weak==0])
sd(cl1.2loneIN$WA[cl1.2loneIN$weak==0])
#sanity check--they are gaps in means but SEs are BIG!#

##no match for weakmr2 so drop that##
cem.loneW<-zelig(formula=cwinit ~ weak + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2loneIN, model="probit", weights=cl1.2loneIN$weights2)
summary(cem.loneW)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.loneW)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2loneIN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2loneIN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.loneW,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.loneW,vcovCL2.2)
ses2.2

##YIELDS ATE OF -0.509 AND SE OF .194 FOR EFFECT OF WEAK WHICH IS SIGNIFICANT BUT IN WRONG DIRECTION##


##SO I THINK THE LARGER AND INSIGNIFICANT SEs ARE RIGHT BUT THEY NEED TO BE ADJUSTED FOR THE NON 1X1 MATCHING!!!


##################################################
##################################################
### CONTROL CASE SHOULD HAVE ONLY 2 BOXES IN IT ##
##################################################
##################################################

##  WEAK HAS AN EFFECT  --> 2 BOXES IN CONTROL AND 2 TREATMENTS##

head(cl1.2)

##Generate treatment and control boxes (agg 5 boxes for control)
cl1.2$treat<-cl1.2$weak
cl1.2$treat[cl1.2$weak==0 & cl1.2$rival==0]<-0
cl1.2$treat[cl1.2$weak==0 & cl1.2$rival==1]<-0
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==0]<-1
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==1]<-2

table(cl1.2$treat)			
head(cl1.2)

##MATCH:: TRY ALL THREE BOXES IN EACH STRATA##
cem.mat3<-cem(treatment="treat", data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival"))
cem.mat3

##LEFT WITH 960C, 66T1, 190T2
960+66+190

head(cl1.2)
cl1.2$IN <- cem.mat3$matched
weights3<-cem.mat3$w
cl1.2<-cbind(cl1.2,weights3)
dim(cl1.2)

cl1.2$WA<-cl1.2$weights3*cl1.2$cwinit
mean(cl1.2$WA[cl1.2$treat==2 & cl1.2$weights3>0])
length(cl1.2$WA[cl1.2$treat==0 & cl1.2$weights3>0])

cl1.2IN <- cl1.2[cl1.2$IN==TRUE,]
dim(cl1.2IN)

mean(cl1.2IN$WA[cl1.2IN$treat==0])
sd(cl1.2IN$WA[cl1.2IN$treat==0])
#sanity check--they are gaps in means but SEs are BIG!#

head(cl1.2IN)

##pol_rel and weakmr2 are all the same in matched data so drop from regression#
cem.W<-zelig(formula=cwinit ~ weak + rival + treat + relcap + jdem + pceyrs + pceyrs2 + pceyrs3, data=cl1.2IN, model="probit", weights=cl1.2IN$weights3)
summary(cem.W)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.W)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.W,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.W,vcovCL2.2)
ses2.2


## SO HERE WE HAVE IT: INTERACTION TERM IS SIGNIFICANT (90% CI) AND LOWER ORDER TERMS ARE ALSO. INTERACTION TERM COEF IS 0.595 (SE 0.3243) BUT SEs ARE WRONG BC OF MATCHING WITH REPLACEMENT ##

##CONFLICTING AGAIN WHERE COARSENED TREATS GET SIGNIFICANCE ON BOTH DICHOT AND TRICHOT WHILE UNCOARSENED ONES DONT AND EVEN WORSE IS THE DICHOT AND TRICHOT SIGNIFICANT RESULTS ARE POINTING IN DIFFERENT DIRECTIONS ON THE INTERACTION TERMS##


#####################################################













#####################################################
#####################################################
## Imai and Van Dyck 2004 Propensity Score Subclassification Strategy: Estimate linear predictor of treatment level, subclassify into a 3x3 matrix, run the probit parametric model in each box to find causal effects, then weight each one by number of observations and get the overall ATE (See Rich's work on this as well)
#####################################################
#####################################################

#Generate Dataset of Interest#
head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)
head(cl1.2)


#LINEAR propensity score for assignment of weakness#
#Don't fit on rival or the interaction term #
w.prop<-zelig(weakmr_decline9 ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2, model="ls")$fitted.values
length(w.prop)
head(w.prop)

#Chop propensity score for weakness into thirds by adding indicator variable for which 3rd it falls into#
sort.w.prop<-sort(w.prop)
cut<-length(sort.w.prop)/3
cut1<-sort.w.prop[cut]
cut1
cut2<-sort.w.prop[2*cut]
cut2
##Will cut for all those <= the cuts##

#Add prop scores to original dataset to keep ordering right for adding prop scores when testing for rivals as well#
w.cut<-w.prop
w.prop<-cbind(w.prop,w.cut)
w.prop[,2][w.prop[,1] <= cut1]<-1
w.prop[,2][w.prop[,1] > cut1 & w.prop[,1] <= cut2]<-2
w.prop[,2][w.prop[,1] > cut2]<-3
w.prop<-as.data.frame(w.prop)
hist(w.prop$w.cut)

##Ok, weak is set, now propensity score for rival## 
r.prop<-zelig(rivalry_count1 ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2, model="ls")$fitted.values
length(r.prop)
head(r.prop)

#Chop propensity score for rival into thirds by adding indicator variable for which 3rd it falls into#
sort.r.prop<-sort(r.prop)
cut<-length(sort.r.prop)/3
cut1<-sort.r.prop[cut]
cut1
cut2<-sort.r.prop[2*cut]
cut2
##Will cut for all those <= the cuts##

#Add prop scores to original dataset to keep ordering right for adding prop scores when testing for rivals as well#
r.cut<-r.prop
r.prop<-cbind(r.prop,r.cut)
r.prop[,2][r.prop[,1] <= cut1]<-1
r.prop[,2][r.prop[,1] > cut1 & r.prop[,1] <= cut2]<-2
r.prop[,2][r.prop[,1] > cut2]<-3
r.prop<-as.data.frame(r.prop)
head(r.prop)
hist(r.prop$r.cut)

## Now lets combine the two cut datasets and place each obs into 1-9 boxes moving across row then down (ie 123 in top row and 789 in bottom row with weak as columns and rivals as rows)##

prop<-cbind(w.prop,r.prop)
head(prop)

prop$box<-prop$r.cut
prop$box[prop$w.cut==1 & prop$r.cut==1]<-1
prop$box[prop$w.cut==2 & prop$r.cut==1]<-2
prop$box[prop$w.cut==3 & prop$r.cut==1]<-3
prop$box[prop$w.cut==1 & prop$r.cut==2]<-4
prop$box[prop$w.cut==2 & prop$r.cut==2]<-5
prop$box[prop$w.cut==3 & prop$r.cut==2]<-6
prop$box[prop$w.cut==1 & prop$r.cut==3]<-7
prop$box[prop$w.cut==2 & prop$r.cut==3]<-8
prop$box[prop$w.cut==3 & prop$r.cut==3]<-9

hist(prop$box)

## Now attach these boxes to the cl1.2 dataset and run the parametric probit model in each of 9 boxes to get ATE and you need robust SEs##

box<-prop$box
cl1.2<-cbind(cl1.2,box)
head(cl1.2)

g1<-cl1.2[cl1.2$box==1,]
g2<-cl1.2[cl1.2$box==2,]
g3<-cl1.2[cl1.2$box==3,]
g4<-cl1.2[cl1.2$box==4,]
g5<-cl1.2[cl1.2$box==5,]
g6<-cl1.2[cl1.2$box==6,]
g7<-cl1.2[cl1.2$box==7,]
g8<-cl1.2[cl1.2$box==8,]
g9<-cl1.2[cl1.2$box==9,]


ate1 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = g1)
ate1<-ate1$coefficients[2:4]
w1<-c(nrow(g1),nrow(g1),nrow(g1))
n<-c(nrow(prop),nrow(prop),nrow(prop))
ate1<-cbind(ate1, w1, n)
ate1<-as.data.frame(ate1)
ate1$wate1<-ate1$ate1*(ate1$w1/ate1$n)
ate1

ate2 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = g2)
ate2<-ate2$coefficients[2:4]
w2<-c(nrow(g2),nrow(g2),nrow(g2))
n<-c(nrow(prop),nrow(prop),nrow(prop))
ate2<-cbind(ate2, w2, n)
ate2<-as.data.frame(ate2)
ate2$wate2<-ate2$ate2 * (ate2$w2/ate2$n)
ate2

ate3 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = g3)
ate3<-ate3$coefficients[2:4]
w3<-c(nrow(g3),nrow(g3),nrow(g3))
n<-c(nrow(prop),nrow(prop),nrow(prop))
ate3<-cbind(ate3, w3, n)
ate3<-as.data.frame(ate3)
ate3$wate3<-ate3$ate3*(ate3$w3/ate3$n)
ate3

ate4 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = g4)
ate4<-ate4$coefficients[2:4]
w4<-c(nrow(g4),nrow(g4),nrow(g4))
n<-c(nrow(prop),nrow(prop),nrow(prop))
ate4<-cbind(ate4, w4, n)
ate4<-as.data.frame(ate4)
ate4$wate4<-ate4$ate4*(ate4$w4/ate4$n)
ate4

ate5 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = g5)
ate5<-ate5$coefficients[2:4]
w5<-c(nrow(g5),nrow(g5),nrow(g5))
n<-c(nrow(prop),nrow(prop),nrow(prop))
ate5<-cbind(ate5, w5, n)
ate5<-as.data.frame(ate5)
head(ate5)
ate5$wate5<-ate5$ate5 * (ate5$w5/ate5$n)
ate5

ate6 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = g6)
ate6<-ate6$coefficients[2:4]
w6<-c(nrow(g6),nrow(g6),nrow(g6))
n<-c(nrow(prop),nrow(prop),nrow(prop))
ate6<-cbind(ate6, w6, n)
ate6<-as.data.frame(ate6)
ate6$wate6<-ate6$ate6*(ate6$w6/ate6$n)
ate6

ate7 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = g7)
ate7<-ate7$coefficients[2:4]
w7<-c(nrow(g7),nrow(g7),nrow(g7))
n<-c(nrow(prop),nrow(prop),nrow(prop))
ate7<-cbind(ate7, w7, n)
ate7<-as.data.frame(ate7)
ate7$wate7<-ate7$ate7*(ate7$w7/ate7$n)
ate7

ate8 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = g8)
ate8<-ate8$coefficients[2:4]
w8<-c(nrow(g8),nrow(g8),nrow(g8))
n<-c(nrow(prop),nrow(prop),nrow(prop))
ate8<-cbind(ate8, w8, n)
ate8<-as.data.frame(ate8)
ate8$wate8<-ate8$ate8*(ate8$w8/ate8$n)
ate8

ate9 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = g9)
ate9<-ate9$coefficients[2:4]
w9<-c(nrow(g9),nrow(g9),nrow(g9))
n<-c(nrow(prop),nrow(prop),nrow(prop))
ate9<-cbind(ate9, w9, n)
ate9<-as.data.frame(ate9)
ate9$wate9<-ate9$ate9*(ate9$w9/ate9$n)
ate9

## Now aggregate the weighted average treatment effect which in my coding is simply summing the "wate" for each treatment variable ##


ate.rival<-sum(ate1[1,4],ate2[1,4],ate3[1,4],ate4[1,4],ate5[1,4],ate6[1,4],ate7[1,4],ate8[1,4],ate9[1,4])
ate.rival

ate.weak<-sum(ate1[2,4],ate2[2,4],ate3[2,4],ate4[2,4],ate5[2,4],ate6[2,4],ate7[2,4],ate8[2,4],ate9[2,4])
ate.weak

ate.inter<-sum(ate1[3,4],ate2[3,4],ate3[3,4],ate4[3,4],ate5[3,4],ate6[3,4],ate7[3,4],ate8[3,4],ate9[3,4])
ate.inter


### SO HERE IT IS: Weighted Avg Treatment Effect rival(-0.0082) and weak (-0.1521) and interacted ( 0.0524). BUT WE DON'T HAVE SEs###########





############################################################
############################################################
## OK, NOW LET'S TRY ACTING AS IF THE INTERACTION TERM IS ACTUAL THE TREATMENT ITSELF (WHICH IS WHAT THE AUTHORS SEEM TO SUGGEST) AND SIMPLY COARSEN ON DIFFERING VALUES OF NRIV. THE QUESTION WE HAVE HERE CONCERNS WHETHER TO INCLUDE THE LOWER ORDER TERMS OR DO WE SIMPLY THINK OF NRIV AS SOME NEW AGGREGATED STANDALONE TERM AND IF WE DO INCLUDE IT DO WE NEED BALANCE ON IT OR LET IT FLOAT##
############################################################
############################################################

#Generate Dataset of Interest#
head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)
head(cl1.2)

#Run Naive Probit Model 1.2 without Lower Order Terms#
int <- zelig(cwinit ~ nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = cl1.2)


mat2.2<-estfun(int)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(int,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(int,vcovCL2.2)
ses2.2

## Naive model without lower order terms suggests interaction itself isn't meaningful w/ coef of 0.0186 and SE of 0.0177 ##

#####################################

##Now coarsen interaction term (NRIV) and estimate WITHOUT lower order terms ##
head(cl1.2)
##4 BUCKETS##
hist(cl1.2$nriv_weakdec9mr)
cl1.2$int<-cl1.2$dyad
cl1.2$int[cl1.2$nriv==0]<-1
cl1.2$int[cl1.2$nriv>0 & cl1.2$nriv <=2]<-2
cl1.2$int[cl1.2$nriv>2 & cl1.2$nriv <=3.5]<-3
cl1.2$int[cl1.2$nriv>3.5]<-4
hist(cl1.2$int)

##MATCH:: TRY ALL FOUR BOXES IN EACH STRATA##
cem.int<-cem(treatment="int", data=cl1.2, drop=c("cwinit","nriv_weakdec9mr","rivalry_count1","weakmr_decline9","dyad"))
cem.int

##LEFT WITH 1274C, 136T1, 133T2, 77T3
1274+136+133+77

head(cl1.2)
cl1.2$IN <- cem.int$matched
weights3<-cem.int$w
cl1.2<-cbind(cl1.2,weights3)
dim(cl1.2)

cl1.2$WA<-cl1.2$weights3*cl1.2$cwinit
mean(cl1.2$WA[cl1.2$int==4 & cl1.2$weights3>0])
length(cl1.2$WA[cl1.2$int==4 & cl1.2$weights3>0])

cl1.2IN <- cl1.2[cl1.2$IN==TRUE,]
dim(cl1.2IN)

mean(cl1.2IN$WA[cl1.2IN$int==4])
sd(cl1.2IN$WA[cl1.2IN$int==4])
#sanity check--they are gaps in means but SEs are BIG!#

head(cl1.2IN)

##weakmr2 are all the same in matched data so drop from regression#
int.W<-zelig(formula=cwinit ~ nriv_weakdec9mr + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2IN, model="probit", weights=cl1.2IN$weights3)
summary(int.W)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(int.W)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(int.W,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(int.W,vcovCL2.2)
ses2.2

## FIND IN MODEL W/ 4 BUCKETS FOR NRIV WHILE DROPPING LOWER ORDER TERMS AND UNCOARSENING TREATMENT IN FINAL ANALYSIS WE SEE NO SIGNIFICANT EFFECT WITH COEF (0.069) AND SE (0.055)



## NEXT STEP IS RERUN WITH LOWER ORDER TERMS IN UNMATCHED AND THEN AGAIN WITH THEM MATCHED ##



###################################
## Now run model w/ interaction still coarsened into 4 buckets while including lower order terms but leaving them unmatched ##

head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)
head(cl1.2)

hist(cl1.2$nriv_weakdec9mr)
cl1.2$int<-cl1.2$dyad
cl1.2$int[cl1.2$nriv==0]<-1
cl1.2$int[cl1.2$nriv>0 & cl1.2$nriv <=2]<-2
cl1.2$int[cl1.2$nriv>2 & cl1.2$nriv <=3.5]<-3
cl1.2$int[cl1.2$nriv>3.5]<-4
hist(cl1.2$int)

cem.int<-cem(treatment="int", data=cl1.2, drop=c("cwinit","nriv_weakdec9mr","rivalry_count1","weakmr_decline9","dyad"))
cem.int

##LEFT WITH 1274C, 136T1, 133T2, 77T3
1274+136+133+77

head(cl1.2)
cl1.2$IN <- cem.int$matched
weights3<-cem.int$w
cl1.2<-cbind(cl1.2,weights3)
dim(cl1.2)

cl1.2$WA<-cl1.2$weights3*cl1.2$cwinit
mean(cl1.2$WA[cl1.2$int==4 & cl1.2$weights3>0])
length(cl1.2$WA[cl1.2$int==4 & cl1.2$weights3>0])

cl1.2IN <- cl1.2[cl1.2$IN==TRUE,]
dim(cl1.2IN)

mean(cl1.2IN$WA[cl1.2IN$int==4])
sd(cl1.2IN$WA[cl1.2IN$int==4])
#sanity check--they are gaps in means but SEs are BIG!#

head(cl1.2IN)

##weakmr2 are all the same in matched data so drop from regression#
int.W<-zelig(formula=cwinit ~ nriv_weakdec9mr + rivalry_count1 + weakmr_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2IN, model="probit", weights=cl1.2IN$weights3)
summary(int.W)
##getting error message but is generating something##


##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(int.W)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(int.W,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(int.W,vcovCL2.2)
ses2.2

## Now get significance on both the interaction term and weakness term: interaction (0.413 and SE .135) ##

#################################################

## Now run model w/ interaction still coarsened into 4 buckets while including lower order terms and matching on them ##

head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)
head(cl1.2)

hist(cl1.2$nriv_weakdec9mr)
cl1.2$int<-cl1.2$dyad
cl1.2$int[cl1.2$nriv==0]<-1
cl1.2$int[cl1.2$nriv>0 & cl1.2$nriv <=2]<-2
cl1.2$int[cl1.2$nriv>2 & cl1.2$nriv <=3.5]<-3
cl1.2$int[cl1.2$nriv>3.5]<-4
hist(cl1.2$int)

cem.int<-cem(treatment="int", data=cl1.2, drop=c("cwinit","nriv_weakdec9mr","dyad"))
cem.int


## CANNOT BE DONE BECAUSE MATCHES WON'T EXIST IF BALANCING LOWER ORDER TERMS WHILE LOOKING IN DIFFERENT BUCKETS OF INTERACTION TERM ##


#####################################################
#####################################################
##       WEAK ON WEAK THESIS   ######################
#####################################################
#####################################################

## Now let's check out WEAK ON WEAK THESIS ##

# OK, so we already know that a weak target on its own makes dispute initiation more likely from Model 1.1 w/ no interaction#

## (1) Naive regression w/ interaction btw weak and weak ##

head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)
head(cl1.2)

#Dichotomize having a weak target and being weak yourself#
cl1.2$wt<-cl1.2$weakmr2_decline9
cl1.2$wt[cl1.2$weakmr2_decline9 == 0]<-0
cl1.2$wt[cl1.2$weakmr2_decline9 > 0]<-1
cl1.2$ww<-cl1.2$weakmr_decline9*cl1.2$wt

##Invert weakness yourself: treatment is now, you are strong and the target is weak#
cl1.2$wy<-cl1.2$weakmr_decline9
cl1.2$wy[cl1.2$weakmr_decline9 == 0]<-1
cl1.2$wy[cl1.2$weakmr_decline9 > 0]<-0
cl1.2$ww2<-cl1.2$wy*cl1.2$wt

##Invert weakness but make continuous ##
cl1.2$strong<-1 - cl1.2$weakmr_decline9
head(cl1.2)



##Descriptive##
mean(cl1.2$cwinit[cl1.2$wy == 0 & cl1.2$wt ==0])
mean(cl1.2$cwinit[cl1.2$wy == 0 & cl1.2$wt ==1])
mean(cl1.2$cwinit[cl1.2$wy == 1 & cl1.2$wt ==0])
mean(cl1.2$cwinit[cl1.2$wy == 1 & cl1.2$wt ==1])
mean(cl1.2$cwinit[cl1.2$wy == 0])
mean(cl1.2$cwinit[cl1.2$wy == 1])
mean(cl1.2$cwinit[cl1.2$wt == 0])
mean(cl1.2$cwinit[cl1.2$wt == 1])
## LOOKS LIKE THE MOST OBVIOUS STORY IS THE ONE GOING ON, STATES WITH STRONG REPS ARE INITIATING DISPUTES WITH STATES WITH WEAK REPS IF THEY HAVE BEEN IN A DISPUTE RECENTLY ##
rivalry_count1 +  +   + pceyrs + pceyrs2 + pceyrs3

ww<-zelig(formula=cwinit ~ strong + wt + strong*wt + relcap + jdem + pol_rel, data=cl1.2, model="probit")
summary(ww)


x.low<-setx(ww, wt=0)
x.high<-setx(ww, wt=1)
s.out<- sim(ww,x=x.low,x1=x.high)
summary(s.out)
plot(s.out)



mat2.2<-estfun(ww)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(ww,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(ww,vcovCL2.2)
ses2.2


#####################################
## BASIC THESIS: WEAK REP GETS YOU ##
## ATTACKED W/ COARSENING & #########
## NO INTERACTION TERMS   ###########
#####################################


# Generate dataset of interest Model 1.1 #
head(cla)
cl1.1<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.1)
head(cl1.1)

#Resetting rivalry count 1 and the interaction term
cl1.1$rivalry_count1_orig<-cl1.1$rivalry_count1
cl1.1$rivalry_count1<-cl1.1$rivalry_count1_orig-1

cl1.1$nriv_weakdec9mr_orig<-cl1.1$nriv_weakdec9mr
cl1.1$nriv_weakdec9mr<-cl1.1$rivalry_count1 * cl1.1$weakmr_decline9

#Drop original terms here
dim(cl1.1)
cl1.1<-cl1.1[,-13]
cl1.1<-cl1.1[,-12]
head(cl1.1)



## Coarsen weak target into 2 buckets ##

cl1.1$wktar<-cl1.1$weakmr2_decline9
cl1.1$wktar[cl1.1$weakmr2_decline9==0]<-1
cl1.1$wktar[cl1.1$weakmr2_decline9>0 ]<-2

sort(unique(cl1.1$weakmr2_decline9))
hist(cl1.1$wktar)

##Create Strong Challenger Variable
cl1.1$strong<-1-cl1.1$weakmr_decline9
cl1.1$strong[cl1.2$strong<1]<-0
head(cl1.1)
hist(cl1.1$strong)
hist(cl1.1$weakmr2_decline9)

##Creat treat var with strong and weak target##
cl1.1$treat<-cl1.1$strong
cl1.1$treat[cl1.1$wktar==1 & cl1.1$strong==0]<-0
cl1.1$treat[cl1.1$wktar==1 & cl1.1$strong==1]<-0
cl1.1$treat[cl1.1$wktar==2 & cl1.1$strong==0]<-0
cl1.1$treat[cl1.1$wktar==2 & cl1.1$strong==1]<-1

## Run CEM ##
head(cl1.1)
cem.int<-cem(treatment="wktar", data=cl1.1, drop=c("cwinit","dyad","weakmr2_decline9","strong","wktar"))
cem.int


##LEFT WITH 
2668+1075

head(cl1.1)
cl1.1$IN <- cem.int$matched
weights<-cem.int$w
cl1.1<-cbind(cl1.1,weights)
dim(cl1.1)

cl1.1$WA<-cl1.1$weights
cl1.1$WA<-cl1.1$weights*cl1.1$cwinit

#head(cl1.1)
#mean(cl1.2$WA[cl1.2$int==4 & cl1.2$weights3>0])
#length(cl1.2$WA[cl1.2$int==4 & cl1.2$weights3>0])

cl1.1IN <- cl1.1[cl1.1$IN==TRUE,]
dim(cl1.1IN)

mean(cl1.1IN$WA[cl1.1IN$wktar==3])
sd(cl1.1IN$WA[cl1.1IN$int==4])
#sanity check--they are gaps in means 

summary(cl1.1IN$treat)

##weakmr are all the same in matched data so drop from regression#
int<-zelig(formula=cwinit ~ weakmr2_decline9 + weakmr_decline9 + rivalry_count1 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.1IN, model="probit", weights=cl1.1IN$weights)
summary(int)


mat2.2<-estfun(int)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.1IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.1IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.4<-df2.2*sandwich(int,meat=crossprod(u2.2)/N2.2)
(vcovCL2.4)
ses2.4<-coeftest(int,vcovCL2.4)
ses2.4



#############################################
## TESTING WEAK REPUTATION OF TARGET:: WITH INTERACTION ##############################################
##############
cl.wk2<-na.omit(cla[c("cwinit","rivalry_count2","weakmr_decline9","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])

#Reset rivalry
#Resetting rivalry count 2 
cl.wk2$rivalry_count2_orig<-cl.wk2$rivalry_count2
cl.wk2$rivalry_count2<-cl.wk2$rivalry_count2_orig-1
#Delete excess rival count column
cl.wk2<-cl.wk2[,-12]

#Run naive probit
wk2<-zelig(cwinit~rivalry_count2+weakmr_decline9+weakmr2_decline9+relcap+jdem+pol_rel+pceyrs+pceyrs2+pceyrs3 + rivalry_count2*weakmr2_decline9, data=cl.wk1, model="probit")
summary(wk2)

#Run sandwich
mat2.2<-estfun(wk2)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl.wk2$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl.wk2$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcov.wk2<-df2.2*sandwich(wk2,meat=crossprod(u2.2)/N2.2)
(vcov.wk2)
ses.wk2<-coeftest(wk2,vcov.wk2)
ses.wk2

table.wk2<-ses.wk2[,1:2]
xtable(table.wk2,digits=3,caption="naive weak rep (with interaction)")
## FIND INTERACTION ISN'T SIGNIFICANT BUT WEAK REP LOWER ORDER IS SIGNIFICANT TO GETTING TARGETED (0.368 and SE 0.148)!

##### NOW MATCH FOR 1.2 WEAK TARGET #####################

##NOW CEM INTERACTION INTO TWO BUCKETS THEN RE-RUN#
head(cl.wk2)
#generate weak2 rep dummy#
cl.wk2$int<-cl.wk2$weakmr2_decline9
cl.wk2$int[cl.wk2$weakmr2_decline9==0 | cl.wk2$rivalry_count2==0]<-0
cl.wk2$int[cl.wk2$weakmr2_decline9>0 & cl.wk2$rivalry_count2>0]<-1

cem.wk2<-cem(treatment="int", data=cl.wk2, drop=c("cwinit","weakmr2_decline9","dyad","int","rivalry_count2"))
cem.wk2
## Keep 4813 Control and 1202 Treated with L1 of 0.593 ##

#Add variable noting whether that observation was matched and what weight to put on it
cl.wk2$IN <- cem.wk2$matched
weight.wk2<-cem.wk2$w
cl.wk2<-cbind(cl.wk2,weight.wk2)
dim(cl.wk2[cl.wk2$IN==TRUE,])

#Create dataset with only matched observation
cl.wk2IN <- cl.wk2[cl.wk2$IN==TRUE,]
dim(cl.wk2IN)

#Run parametric model on the matched dataset with proper weights
cem.wk2.M<-zelig(formula=cwinit ~ rivalry_count2 + weakmr_decline9 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3 + rivalry_count2*weakmr2_decline9, data=cl.wk2IN, model="probit", weights=cl.wk2IN$weight.wk2)
##WHAT IS THIS ERROR MESSAGE HERE??
summary(cem.wk2.M)

#Get clustered standard errors
mat2.2<-estfun(cem.wk2.M)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl.wk2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl.wk2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcov.wk2.M<-df2.2*sandwich(cem.wk2.M,meat=crossprod(u2.2)/N2.2)
(vcov.wk2.M)
ses.wk2.M<-coeftest(cem.wk2.M,vcov.wk2.M)
ses.wk2.M

table.wk2.M<-ses.wk2.M[,1:2]
xtable(table.wk2.M,digits=3,caption="weak2 rep matched has an effect BUT INT DOESN'T")
##NO EFFECT FOR INTERACTION POST-MATCHING BUT SOME FOR WEAK2 REP

























#######################################
## 1st step: test effect of ###########
## increasing rivals while ############
## holding rep at strong  #############
#######################################

cl1.2strong<-cl1.2[cl1.2$weak==0,]
dim(cl1.2strong)

### Naive probit regression ############
########################################
## IN ALL OF THE BELOW REGRESSIONS I HAVE TO PLUCK OUT VARIABLES THAT HAVE NOT VARIATION (WHICH EVER OF THE TREATMENTS IS FIXED AT ZERO PLUS THOSE THAT POST-MATCHING HAVE ONLY ONE VALUE LEFT) OTHERWISE THE SANDWICH FUNCTION CHOKES ####
#######################################



naive.strong<-zelig(formula=cwinit ~ rivalry_count1 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2strong, model="probit")
summary(naive.strong)

mat2.2<-estfun(naive.strong)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2strong$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2strong$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(naive.strong,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(naive.strong,vcovCL2.2)
ses2.2

length(N2.2)
dim(u2.2)

## Naive probit varying level of rivals with a 0 rep yields insignficant effect for rivalry count (-0.014) and std (0.018) ##

#########################################

## Match each level of rivals to control
## level of rivals then estimate separate 
## ATE for each one ####################

##trying to do multi treatments at once##
##SWITCHED TO LOGIT BC PROBIT NOT IMPLEMENTED IN CEM ATT YET##


unique(cl1.2strong$rival)

cem.mat<-cem(treatment="rival", data=cl1.2strong, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.mat$tab
summary(cem.mat)

cem.strong<-att(obj=cem.mat,formula=cwinit ~ rivalry_count1 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2strong,model="logit")
summary(cem.strong)
names(cem.strong)
cem.strong$att.model
class(cem.strong$att.model)
cem.strong<-as.glm(cem.strong)



###CAN'T GET CLUSTERED STD ERRORS BC ATT PUMPS OUT A NON-LM OBJECT###

mat2.2<-estfun(cem.strong)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2strong$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2strong$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.strong,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.strong,vcovCL2.2)
ses2.2


##GET 1143, 1853, 1467 MATCHED UNITS##
##RIVALRY COUNT ATT = .044 BUT NO READ ON SEs

##let's try weighting##
cem.strongW<-zelig(formula=cwinit ~ rivalry_count1 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2strong, model="probit", weights=cem.mat$w)
summary(cem.strongW)
##getting error message but is generating something##

mat2.2<-estfun(cem.strongW)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2strong$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2strong$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.strongW,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.strongW,vcovCL2.2)
ses2.2

##YIELDS ATE OF .0216 AND SE OF .0143 FOR EFFECT OF RIVALS##




















####################
####################

##trying to do multi treatments  ATT at once##
##SWITCHED TO LOGIT BC PROBIT NOT IMPLEMENTED IN CEM ATT YET##
cl1.2lone<-cl1.2[cl1.2$rival==0,]
dim(cl1.2lone)
head(cl1.2lone)

unique(cl1.2lone$weak)

cem.mat2<-cem(treatment="weak", data=cl1.2lone, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.mat2
##PINNING DOWN THE ACTUAL DATASET POST-MATCHING BUT THE WEIGHTS AREN'T ACCOUNTED FOR##
cem.mat2$matched[cem.mat2$matched==TRUE] <-1
cem.matched2<-cem.mat2$matched
cl1.2lone<-cbind(cl1.2lone,cem.matched2)
cl1.2loneM<-cl1.2lone[cl1.2lone$cem.matched2==1,]
dim(cl1.2loneM)

cem.mat2$w[cem.mat2$w==0]<-NA
cem.weight2<-na.omit(cem.mat2$w)
cl1.2loneMW<-cbind(cl1.2loneM,cem.weight2)
dim(cl1.2loneMW)

mean(cl1.2$cwinit[cl1.2$weak==2])

cem.loneMW<-zelig(formula=cwinit ~ weakmr_decline9 + relcap + jdem + pceyrs + pceyrs2 + pceyrs3, data=cl1.2loneMW, model="probit", weights=cl1.2loneMW$cem.weight2)
summary(cem.loneMW)

mat2.2<-estfun(cem.loneMW)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2loneMW$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2loneMW$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.loneMW,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.loneMW,vcovCL2.2)
ses2.2



########
cem.lone<-att(obj=cem.mat2,formula=cwinit ~ weakmr_decline9 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2lone,model="logit")
summary(cem.lone)
names(cem.lone)
cem.lone$att.model
class(cem.lone$att.model)
cem.lone<-as.glm(cem.lone)




###CAN'T GET CLUSTERED STD ERRORS BC ATT PUMPS OUT A NON-LM OBJECT###

mat2.2<-estfun(cem.strong)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2strong$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2strong$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.strong,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.strong,vcovCL2.2)
ses2.2

##GET 116, 32, 20 MATCHED UNITS##
##RIVALRY COUNT ATT = -0.677 BUT NO READ ON SEs

weights2<-cem.mat2$w
cl1.2lone<-cbind(cl1.2lone,weights2)
dim(cl1.2lone)

cl1.2lone$WA<-cl1.2lone$weights2*cl1.2lone$cwinit
mean(cl1.2lone$WA[cl1.2lone$weak==2 & cl1.2lone$weights2>0])
length(cl1.2lone$WA[cl1.2lone$weak==0 & cl1.2lone$weights2>0])
##OK SANITY CHECK WORKS--SIG EFFECT GOING FROM 0 to 1###

##let's try weighting##
cem.loneW<-zelig(formula=cwinit ~ weakmr_decline9 + relcap + jdem + pceyrs + pceyrs2 + pceyrs3, data=cl1.2lone, model="probit", weights=cl1.2lone$weights2)
summary(cem.loneW)
##getting error message but is generating something##


































####################

## Create dataset with with rival=0 & =1

cl1.2strong1<-cl1.2strong[cl1.2strong$rival<2,]
dim(cl1.2strong1)
unique(cl1.2strong$weak)

strong1<-matchit(formula=rival ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3,data=cl1.2strong1, method="cem")
summary(strong1)
## matches 1633 C to 2511 T1 ##
strong1
strong1.data<-match.data(strong1)
dim(strong1.data)


pre<-summary(strong1)$sum.all
post<-summary(strong1)$sum.match

pre.std<-(pre[,1]-pre[,2])/pre[,3]
post.std<-(post[,1]-post[,2])/post[,3]
table<-cbind(pre.std,post.std)
table


plot(x = pre.std, y = 1:nrow(pre),
     pch = 19, col = "burlywood2", xlab = "Standardized Imbalance: CEM Strong 1", axes = FALSE, 
     xlim = c(-3,3), ylab = "")
points(x = post.std, y = 1:nrow(post),
     pch = 19, col = "cadetblue2")
abline(v = 0, lty = "dashed")
axis(1)
axis(2, at = 1:8, labels = rownames(pre), las = 1, cex.axis = .8)
legend("topright", legend = c("Pre-match","Post-match"), pch = c(19,19),
   col = c("burlywood2","cadetblue2"), cex = .8)

##NEED TO DROP THE ORIG VARIABLES FROM THE IMBALANCE TESTS###
pre.imb<-imbalance(group=cl1.2strong1$rival, data=cl1.2strong1, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival"))

post.imb<-imbalance(group=strong1.data$rival, data=strong1.data, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
pre.imb
post.imb
##REALLY BAD L1 despite great univariate balance ##

cem.strong1<-zelig(formula=cwinit ~ rivalry_count1 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=strong1.data, model="probit")
summary(cem.strong1)
cem.strong1
unique(strong1.data$weakmr2_decline9)

class(cem.strong1)

mat2.2<-estfun(cem.strong1)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,strong1.data$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(strong1.data$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.strong1,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.strong1,vcovCL2.2)
ses2.2


############# NEED TO FIGURE THIS OUT ##################
cem.match<-cem(treatment="rival", data=cl1.2strong1, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
cem.match
cem.match$n.strata
summary(cem.match)
cem.match$w[1:50]

cem.strong1a<-zelig(formula=cwinit ~ rivalry_count1 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2strong1, model="probit", weights=cem.match$w)
summary(cem.strong1a)


#######################################################


cem.att<-att(obj=cem.match, formula=cwinit ~ rivalry_count1 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2strong1, model="probit")
##probit isn't implemented yet## (yikes)
##what does this att function do that normal one doesn't and why is L1
##different even w/ seemingly same observations (same # at least)
#######################################################


## Matching all strong rival control to rival treat 1 yields no significance for rival treat 1 (0.0016 w/ std 0.0477) ##

#######################################################

## Create dataset with with rival=0 & =2

cl1.2strong2<-cl1.2strong[cl1.2strong$rival==0 | cl1.2strong$rival==2,]
##have to dichotomize treatment##
cl1.2strong2$rival[cl1.2strong2$rival==2]<-1
dim(cl1.2strong2)
head(cl1.2strong2)

strong2<-matchit(formula=rival ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3,data=cl1.2strong2, method="cem")
summary(strong2)
## matches 1383 C to 1684 T1 ##
strong2
strong2.data<-match.data(strong2)
dim(strong2.data)


pre<-summary(strong2)$sum.all
post<-summary(strong2)$sum.match

pre.std<-(pre[,1]-pre[,2])/pre[,3]
post.std<-(post[,1]-post[,2])/post[,3]
table<-cbind(pre.std,post.std)
table


plot(x = pre.std, y = 1:nrow(pre),
     pch = 19, col = "burlywood2", xlab = "Standardized Imbalance: CEM Strong 2", axes = FALSE, 
     xlim = c(-3,3), ylab = "")
points(x = post.std, y = 1:nrow(post),
     pch = 19, col = "cadetblue2")
abline(v = 0, lty = "dashed")
axis(1)
axis(2, at = 1:8, labels = rownames(pre), las = 1, cex.axis = .8)
legend("topright", legend = c("Pre-match","Post-match"), pch = c(19,19),
   col = c("burlywood2","cadetblue2"), cex = .8)

pre.imb<-imbalance(group=cl1.2strong2$rival, data=cl1.2strong2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival"))

post.imb<-imbalance(group=strong2.data$rival, data=strong2.data, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
pre.imb
post.imb
##PRETTY BAD L1 despite great univariate balance ##

cem.strong2<-zelig(formula=cwinit ~ rivalry_count1 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=strong2.data, model="probit")
summary(cem.strong2)

mat2.2<-estfun(cem.strong2)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,strong2.data$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(strong2.data$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.strong2,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.strong2,vcovCL2.2)
ses2.2



## Matching all strong rival control to rival treat 2 yields no significance for rival treat 2 (0.0041 w/ std 0.0225) ##

####################################################
### WITH REP FIXED AT 0 (IE STRONG), MOVING RIVALS #
### HAS NO EFFECT ON ON OUTCOME--> ALL OF THESE ####
### CASES ARE PART OF CONTROL AS RIVAL TREATMENT ###
### IS NOT "ACTIVATED" IF NO REPUATITIONAL #########
### TREATMENT IS PRESENT AS WELL ###################
####################################################


#######################################
## 1st step: test effect of ###########
## increasing weak rep while ##########
## holding rivals at zero #############
#######################################

cl1.2lone<-cl1.2[cl1.2$rival==0,]
dim(cl1.2lone)
head(cl1.2lone)
### Naive probit regression ############

naive.lone<-zelig(formula=cwinit ~ weakmr_decline9 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2lone, model="probit")
summary(naive.lone)

length(cl1.2lone$dyad)
dim(u2.2)
dim(mat2.2)

mat2.2<-estfun(naive.lone)
mat2.2<-na.omit(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2lone$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2lone$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(naive.lone,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(naive.lone,vcovCL2.2)
ses2.2


## Naive probit varying level of weak with a 0 rivals yields a signficant effect for rep (-0.8375) and std (.296) WHICH IS SIGNIFICANT IN THE WRONG DIRECTION ##
########################################################
#########################################

## Match each level of weak to control
## level of weak then estimate separate 
## ATE for each one ####################

## Create dataset with with weak=0 & =1

cl1.2lone1<-cl1.2lone[cl1.2lone$weak<2,]
dim(cl1.2lone1)
head(cl1.2lone1)
unique(cl1.2lone1$weak)

lone1<-matchit(formula=weak ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3,data=cl1.2lone1, method="cem")
summary(lone1)
## matches 323 C to 65 T1 ##
lone1
lone1.data<-match.data(lone1)
dim(lone1.data)
head(lone1.data)

lone1.data$pol_rel

pre<-summary(lone1)$sum.all
post<-summary(lone1)$sum.match

pre.std<-(pre[,1]-pre[,2])/pre[,3]
post.std<-(post[,1]-post[,2])/post[,3]
post.std[2] <- 0
post.std[5] <- 0
table<-cbind(pre.std,post.std)
table


plot(x = pre.std, y = 1:nrow(pre),
     pch = 19, col = "burlywood2", xlab = "Standardized Imbalance: CEM Lone 1", axes = FALSE, 
     xlim = c(-3,3), ylab = "")
points(x = post.std, y = 1:nrow(post),
     pch = 19, col = "cadetblue2")
abline(v = 0, lty = "dashed")
axis(1)
axis(2, at = 1:8, labels = rownames(pre), las = 1, cex.axis = .8)
legend("topright", legend = c("Pre-match","Post-match"), pch = c(19,19),
   col = c("burlywood2","cadetblue2"), cex = .8)

pre.imb<-imbalance(group=cl1.2lone1$weak, data=cl1.2lone1, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival"))

post.imb<-imbalance(group=lone1.data$weak, data=lone1.data, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
pre.imb
post.imb
##BIG IMPROVEMENT ON L1 BUT STILL BAD ##

cem.lone1<-zelig(formula=cwinit ~ weakmr_decline9 + relcap + jdem +  pceyrs + pceyrs2 + pceyrs3, data=lone1.data, model="probit")
summary(cem.lone1)

mat2.2<-estfun(cem.lone1)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,lone1.data$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(lone1.data$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.lone1,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.lone1,vcovCL2.2)
ses2.2

#sanity check
mean(lone1.data$cwinit[lone1.data$weak==0])
mean(lone1.data$cwinit[lone1.data$weak==1])
length(lone1.data$cwinit[lone1.data$weak==0])

## Matching all lone weak control to weak treat 1 yields significance for weak treat 1 (-1.17 w/ std 0.543) SIGNIFICANT BUT IN THE WRONG DIRECTION!!! ##

#######################################################

## Naive probit varying level of weak with a 0 rivals yields a signficant effect for rep (-0.8375) and std (.365) WHICH IS SIGNIFICANT IN THE WRONG DIRECTION ##
########################################################
#########################################

## Match each level of weak to control
## level of weak then estimate separate 
## ATE for each one ####################

## Create dataset with with weak=0 & =2

cl1.2lone2<-cl1.2lone[cl1.2lone$weak==2 | cl1.2lone$weak==0,]
cl1.2lone2$weak[cl1.2lone2$weak==2]<-1
dim(cl1.2lone2)
head(cl1.2lone2)
unique(cl1.2lone1$weak)

lone2<-matchit(formula=weak ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3,data=cl1.2lone2, method="cem")
summary(lone2)
## matches 202 C to 36 T1 ##
lone2
lone2.data<-match.data(lone2)
dim(lone2.data)
head(lone2.data)

unique(lone2.data$weakmr2_decline9)

pre<-summary(lone2)$sum.all
post<-summary(lone2)$sum.match

pre.std<-(pre[,1]-pre[,2])/pre[,3]
post.std<-(post[,1]-post[,2])/post[,3]
post.std[2] <- 0
table<-cbind(pre.std,post.std)
table


plot(x = pre.std, y = 1:nrow(pre),
     pch = 19, col = "burlywood2", xlab = "Standardized Imbalance: CEM Lone 2", axes = FALSE, 
     xlim = c(-3,3), ylab = "")
points(x = post.std, y = 1:nrow(post),
     pch = 19, col = "cadetblue2")
abline(v = 0, lty = "dashed")
axis(1)
axis(2, at = 1:8, labels = rownames(pre), las = 1, cex.axis = .8)
legend("topright", legend = c("Pre-match","Post-match"), pch = c(19,19),
   col = c("burlywood2","cadetblue2"), cex = .8)

pre.imb<-imbalance(group=cl1.2lone2$weak, data=cl1.2lone2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival"))

post.imb<-imbalance(group=lone2.data$weak, data=lone2.data, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival","weights","subclass"))
pre.imb
post.imb
##BIG IMPROVEMENT ON L1 AND NOW DECENT ##

cem.lone2<-zelig(formula=cwinit ~ weakmr_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=lone2.data, model="probit")
summary(cem.lone2)

mat2.2<-estfun(cem.lone2)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,lone2.data$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(lone2.data$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.lone2,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.lone2,vcovCL2.2)
ses2.2

#sanity check
mean(lone2.data$cwinit[lone2.data$weak==1])
mean(lone1.data$cwinit[lone1.data$weak==1])
length(lone1.data$cwinit[lone1.data$weak==0])

## Matching all lone weak control to weak treat 2 yields insignificance for weak treat 2 (-0.59 w/ std 0.64) ##

#######################################################
#######################################################
## SO AFTER ALL TESTS RUNNING ONE 0 TREAT AGAINST THE OTHER AT T1 VS 0 AND T2 VS 0 THE ONLY SIGNIFICANT RESULT IS WITH RIVALS SET TO 0 AND REP GOING FROM STRONG TO WEAK 1 AND THE EFFECT IS IN THE WRONG DIRECTION VS. THE HYPOTHESIS
##  ------>>>> SAYS ODD THINGS ABOUT WHAT IS THE CONTROL!!

#######################################################
#######################################################







#######################################
#######################################
#######################################
#######################################
#######################################
#######################################
a<-(whatif(~ weakmr2_decline9+relcap+jdem+pol_rel+pceyrs+pceyrs2+pceyrs3,data=cl1.2[cl1.2$weakmr_decline9 > 0 & cl1.2$rivalry_count1 > 1,],cfact=cl1.2[cl1.2$weakmr_decline9 > 0 & cl1.2$rivalry_count1 == 1 ,]))
summary(a)
hist(cl1.2$pceyrs3)



#dropping pceyrs2 and 3
a<-(whatif(~ weakmr2_decline9+relcap+jdem+pol_rel+pceyrs,data=cl1.2[cl1.2$weakmr_decline9 > 0 & cl1.2$rivalry_count1 ==1,],cfact=cl1.2[cl1.2$weakmr_decline9 > 0 & cl1.2$rivalry_count1 > 1 ,]))
summary(a)




b<-(whatif(~ weakmr2_decline9+relcap+jdem+pol_rel,data=cl1.3[cl1.3$weakmr_decline9 >= 0.05 & cl1.3$rivalry_count1 >=2,],cfact=cl1.3[cl1.3$weakmr_decline9 == 0 & cl1.3$rivalry_count1==1,]))
summary(b)


#################################
## CASE WITHIN THE DATA #########
#################################
head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad","ccode1","ccode2","year")])
dim(cl1.2)

## So these are cases where a state w/ weak rep initiates a dispute with the US ##
us.back<-cl1.2[cl1.2$ccode2==220,]
dim(us.back)
us.back<-us.back[us.back$cwinit==1,]
dim(us.back)
us.back<-us.back[us.back$weakmr_decline9>0,]
dim(us.back)
us.back

## Incidents ##
#Russia 1948
#Russia 1963
#Russia 1964
#Russia 1966
#Russia 1967
#Cuba 1971

###################################
###################################
###################################
###################################
#### END MESS AROUND TIME #########
###################################
###################################
###################################


##COEFFICIENTS FOR CONVEX HULL ONLY DATA MODEL##
z.out<-zelig(force~rivalry_count1+weakmr_decline9+nriv_weakdec9mr+weakmr2_decline9+relcap+jdem+pol_rel,data=cl.sub,model="probit")
summary(z.out)


##CLUSTERED STANDARD ERRORS FOR CONVEX HULL ONLY DATA MODEL###

mat2.2<-estfun(z.out)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl.sub$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl.sub$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(z.out,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(z.out,vcovCL2.2)
ses2.2

##TABLE 1.4 ORIGINAL######################
head(cla)
cl1.4<-na.omit(cla[c("force","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","dyad")])

cla.pr.1.4<-zelig(force~rivalry_count1+weakmr_decline9+nriv_weakdec9mr+weakmr2_decline9+relcap+jdem+pol_rel, data=cl1.4, model="probit")
summary(cla.pr.1.4)

##TABLE 1.4 ORIGINAL CLUSTERED STANDARD ERRORS#######
mat2.1<-estfun(cla.pr.1.4)
mat2.1<-na.omit(mat2.1)
dim(mat2.1)
N2.1<-nrow(mat2.1)
u2.1<-apply(mat2.1,2,function(x) tapply (x,cl1.4$dyad,sum))
u2.1<-na.omit(u2.1)
clust2.1<-length(unique(cl1.4$dyad))
a2.1<-clust2.1
b2.1<-nrow(mat2.1)
c2.1<-ncol(mat2.1)
df2.1<-((a2.1)/(a2.1-1))*((b2.1-1)/(b2.1-c2.1))
df2.1
vcovCL2.1<-df2.1*sandwich(cla.pr.1.4,meat=crossprod(u2.1)/N2.1)
(vcovCL2.1)
ses1.4<-coeftest(cla.pr.1.4,vcovCL2.1)
ses1.4

xtable(ses1.4[,1:2],digits=3)


###COMBINED COMPARISON TABLE OF ORIGINAL AND CONVEX HULL ONLY DATA ####
hull1.4<-cbind(ses1.4[,1:2],ses2.2[,1:2])
dim(hull1.4)
hull1.4<-rbind(hull1.4[2:8,],hull1.4[1,])

rownames(hull1.4)<-c("Number of Potential Opponents","Challenger's Weak Reputation","Number of Potential Opponents X Challenger's Weak Reputation","Target's Weak Reputation","Relative Capabilities","Joint Democracy","Politically Relevant","Constant")
xtable(hull1.4,digits=3)



























